Thermal Analysis on Active Heat Dissipation Design with Embedded Flow Channels for Flexible Electronic Devices

Heat generation is a major issue in all electronics, as heat reduces product life, reliability, and performance, especially in flexible electronics with low thermal-conductivity polymeric substrates. In this sense, the active heat dissipation design with flow channels holds great promise. Here, a theoretical model, validated by finite element analysis and experiments, based on the method of the separation of variables, is developed to study the thermal behavior of the active heat dissipation design with an embedded flow channel. The influences of temperature and flow velocity of the fluid on heat dissipation performance were systematically investigated. The influence of channel spacing on heat dissipation performance was also studied by finite element analysis. The study shows that performance can be improved by decreasing the fluid temperature or increasing the flow velocity and channel density. These results can help guide the design of active heat dissipation with embedded flow channels to reduce adverse effects due to excessive heating, thus enhancing the performance and longevity of electronic products.


Introduction
In recent years, electronic products have emphasized intelligence and miniaturization, which require higher performance at the integration level. In particular, flexible electronic devices represent a rapidly developing branch, which has enabled many new applications such as curvilinear electronics and bio-integrated electronics [1][2][3][4][5][6] due to their unique advantages of stretchability and flexibility [7][8][9][10]. Furthermore, hybrid flexible systems have been investigated, such as electrophysiology [11], multimodal [12], large-area pressure [13] and other functions [14]; such functions may lead to overheating in the system. Heat generation is a major issue when these electronic devices are in service. According to the Arrhenius relation, the thermal failure rate of flexible electronic devices has an exponential relationship with temperature [15], in which the failure probability doubles for every 10 • C increase in temperature. Additionally, the use of polymer materials with poor thermal conductivity as the substrate creates more challenges regarding heat dissipation in flexible electronic devices.
In order to solve this special challenge, many thermal designs have been reported to improve heat dissipation in flexible electronic devices. Our group focused on passive thermal designs to control and tune the thermal flow in order to achieve adequate heat dissipation. Li et al. [16] reported on an orthotropic substrate design to control the heat flow paths. Shi et al. [17] proposed a phase change material enabled functional soft composite to serve as the thermal protection substrate for wearable electronics. Other groups have developed various new materials to enhance thermal conductivity. Jeong et al. [18] and Bartlett et al. [19] embedded liquid metal inclusions into a polymer to achieve both flexibility and high thermal conductivity. Sun et al. [20] designed a flexible graphene aerogel-based composite phase change film for personal thermal management applications. Although the above designs were all passive, they hold promise for the thermal management of flexible electronic devices. Compared to passive thermal designs, active designs, which have been extensively studied for silicon electronics, could be more efficient in heat dissipation [21][22][23], and are attracting more and more attention. At present, a substrate with embedded flow channels is the most common active heat dissipation design, where fluid flowing in the channel takes away heat. The earliest microchannel for water on a silicon plate, designed by Tuckeman and Pease in the 1980s, was able to achieve 800 W heat dissipation on a 1 cm 2 chip [24]. Perret et al. [25] performed a numerical analysis of heat transfer in rectangular, rhombic and hexagonal flow channels and found that the former had the lowest thermal resistance. Yong et al. [26] added fins to the traditional rectangular channel and found that a smaller fin angle and smaller fin pitch could better improve heat dissipation performance. Husain and Kim [27] studied the influence of the width-depth ratio of the channel on the heat transfer of the system.
Due to the very poor thermal conductivity of polymer substrates, passive heat dissipation methods cannot provide adequate efficiency. Introducing active fluid flow in channels embedded in soft polymer substrates may significantly improve the efficiency of heat dissipation. However, few research efforts have sought to quantitatively investigate the effects of fluid flowing in the channels embedded in a soft polymer substrate. Because of the complexity and difficulty of solving the active heat dissipation problem of fluid-solid coupling, all the studies mentioned above were primarily numerical analyses, lacking theoretical thermal analyses of the structure with embedded flow channels.
This paper aims to develop a theoretical model to investigate the thermal behavior of an active heat dissipation design with embedded flow channels. Finite element analysis and experiments are also carried out to validate the analytical model. The paper is structured as follows. Section 2 describes the theoretical model. Section 3 presents the results and a discussion. Section 4 presents the conclusion.

Theoretical Model
The active heat dissipation design with flow channels can be applied to the substrate for flexible electronic devices. Figure 1a schematically shows a typical soft substrate, which is made of polydimethylsiloxane (PDMS), with an embedded flow channel for active heat dissipation. Since the diameter of the flow channel is usually much smaller than the characteristic length of the substrate (e.g., thickness), axisymmetric treatment of the problem is reasonable. The hollow channel is located at the center of the substrate. The fluid can flow in the channel along the direction indicated by the arrow changing from blue to red. Among the many liquid materials flowing in the channel, water is the most common, due to its nontoxicity, large specific heat capacity and availability. Figure 1b is a schematic diagram of the theoretical model. The origin of coordinates (r, z) is set at the center of the middle surface of the flow channel. R c and R s denote the channel radius and substrate radius, respectively. H is the length of the flow channel. Path 1, Path 2 and Path 3 represent the inlet surface, the fluid-solid contact surface and the symmetric surface, respectively.
The fluid-solid coupling at the interface (r = Rc) can be simplified as a forced convention, while the coefficient hf of heat convention can be obtained by the Sieder-Tate formula as [28], where f c are the Reynolds number representing the convective strength and the Prandtl number representing relationship between temperature and flow boundary layers of the fluid, respectively, and kf, ρf, μf, cf are the thermal conductivity, density, viscosity and specific heat capacity of the fluid, respectively. Since the radius of the flow channel is smaller than that of substrate, the fluid has a negligible influence on the outer boundary of the PDMS. Therefore, the boundary conditions of the PDMS along the radial direction can be expressed as where k and T0 are the thermal conductivity and initial temperature of the substrate, respectively, and Tf is the temperature of fluid.
The fluid temperature in the channel remains almost the same [29], such that the plane at z = 0 can be considered a symmetric plane. Moreover, a natural convection occurs at the upper surface of the PDMS. Thus, the boundary conditions along the axial direction can be expressed as: The steady temperature field T(r, z) of the PDMS substrate satisfies the following axisymmetric heat transfer equation The fluid-solid coupling at the interface (r = R c ) can be simplified as a forced convention, while the coefficient h f of heat convention can be obtained by the Sieder-Tate formula as [28], where are the Reynolds number representing the convective strength and the Prandtl number representing relationship between temperature and flow boundary layers of the fluid, respectively, and k f , ρ f , µ f , c f are the thermal conductivity, density, viscosity and specific heat capacity of the fluid, respectively. Since the radius of the flow channel is smaller than that of substrate, the fluid has a negligible influence on the outer boundary of the PDMS. Therefore, the boundary conditions of the PDMS along the radial direction can be expressed as where k and T 0 are the thermal conductivity and initial temperature of the substrate, respectively, and T f is the temperature of fluid. The fluid temperature in the channel remains almost the same [29], such that the plane at z = 0 can be considered a symmetric plane. Moreover, a natural convection occurs at the upper surface of the PDMS. Thus, the boundary conditions along the axial direction can be expressed as: where h a is the natural convective coefficient and T a is the ambient temperature. We define the temperature rise from ambient temperature as θ = T − T a . The heat transfer Equation (1) and boundary conditions in Equations (3) and (4) Micromachines 2021, 12, 1165 The method of separation of variables was adopted to solve the above heat transfer problem. The temperature rise θ takes the following form: where Y(z) and X(r) satisfy and with γ as the eigenvalue. Equations (7), (9) and (10) yield the solution of θ as: where I 0 and K 0 is the 0-th order modified Bessel function of the first and second kinds, respectively, γ n satisfies kγ sin(γH/2) − h a cos(γH/2) = 0, and coefficients A n and B n satisfy Here Y(γ n , z)dz, I 1 and K 1 are the first-order modified Bessel function of the first and second kinds, respectively. Coefficients A n and B n can be given by: The temperature field T(r, z) can then be obtained as:

Model Verification
We carried out experiments and finite element analysis (FEA) to validate our theoretical model. Figure 2a shows the experimental setup. A flow channel with a radius R c of 1 mm and a length H of 60 mm was embedded in a PDMS substrate with the radius R s of 25 mm. The inlet was connected to an injector containing ice in order to keep the fluid cold. The cold fluid was pumped to the flow channel by an injection pump at a flow rate of 90 mL/min, and then flowed out through the PDMS into a liquid storage tank. The tem-Micromachines 2021, 12, 1165 5 of 10 perature of the PDMS was measured using an infrared thermal imager. Figure 2b shows the infrared thermal imager measurements at the inlet of channel. The fluid temperature at the inlet was 4.3 • C and the environmental temperature T a was 24.6 • C. Figure 2c shows the temperature distribution on the upper inlet surface of the structure. The fluid-solid coupling FEA was performed using the commercial software, FLU-ENT. The thermal conductivity k of PDMS was 0.17 W/(K·m). The parameters of fluid (i.e., water) were set to 0.6 W/(K·m) for the thermal conductivity kf, 998.2 Kg/m 3 for the mass density ρf, 0.001003 Pa·s for the viscosity μf, and 4182 J/(kg·K) for the specific heat capacity cf. The initial temperature T0 was equal to the ambient temperature Ta; both were set to 24.6 °C. The coefficient of thermal convection ha was defined as 20 W/(K·m 2 ) [30].
In order to measure the heat dissipation capacity of the flow channel, we defined a relative temperature drop affected by flow channel as: Figure 3 shows the temperature along Path 1, where the bar representing the experimental error was obtained from three experiments. The theoretical predictions agreed well with the FEA and experiments, which validates the accuracy of our theoretical model. As expected, the temperature increased as the distance to the channel increased, which yielded a decreasing heat dissipation capacity. For example, at r=2 mm, the temperature of PDMS was 12.3 °C, and = 61.3 % Q . When r increased to 6 mm, the temperature increased to 19.7 °C and Q decreased to = 25.6 % Q . These results indicate that the flow channel significantly reduced the temperature of the PDMS. If a relative temperature drop of 25% was set as a criterion for heat dissipation, the influence range of flow channel The fluid-solid coupling FEA was performed using the commercial software, FLUENT. The thermal conductivity k of PDMS was 0.17 W/(K·m). The parameters of fluid (i.e., water) were set to 0.6 W/(K·m) for the thermal conductivity k f , 998.2 Kg/m 3 for the mass density ρ f , 0.001003 Pa·s for the viscosity µ f , and 4182 J/(kg·K) for the specific heat capacity c f . The initial temperature T 0 was equal to the ambient temperature T a ; both were set to 24.6 • C. The coefficient of thermal convection h a was defined as 20 W/(K·m 2 ) [30].
In order to measure the heat dissipation capacity of the flow channel, we defined a relative temperature drop affected by flow channel as: Figure 3 shows the temperature along Path 1, where the bar representing the experimental error was obtained from three experiments. The theoretical predictions agreed well with the FEA and experiments, which validates the accuracy of our theoretical model. As expected, the temperature increased as the distance to the channel increased, which yielded a decreasing heat dissipation capacity. For example, at r = 2 mm, the temperature of PDMS was 12.3 • C, and Q= 61.3%. When r increased to 6 mm, the temperature increased to 19.7 • C and Q decreased to Q= 25.6%. These results indicate that the flow channel significantly reduced the temperature of the PDMS. If a relative temperature drop of 25% was set as a criterion for heat dissipation, the influence range of flow channel on temperature was six times its own radius R c .

Influences of Fluid Temperature and Flow Velocity
In order to investigate the influence of fluid temperature on the thermal behavior of the active heat dissipation design, theoretical predictions and FEA were carried out with Tf of 5/10/15 °C. The initial and ambient temperature were 25 °C and the fluid velocity was 100 mm/s. Figure 4a shows the temperature along Path 2. Again, there was good agreement between the theoretical prediction and FEA and the maximum temperature difference of 0.56 °C. As expected from FEA, the temperature of the PDMS decreased slightly along the flowing direction, since the heat was taken away by the flowing fluid. This slight temperature change validated the theoretical treatment of the influence of flowing fluid as a forced convection with a constant convective coefficient, although the local convective coefficient decreased along the flowing direction. Figure 4b compares the temperature along Path 3 according to the theoretical prediction and FEA with various fluid temperatures of 5/10/15 °C. It can be observed that the theoretical prediction agreed well with FEA. As expected, lower fluid temperature yielded better heat dissipation performance. At the location of r = 10 mm, the temperatures of PDMS were 19.6 °C, 21 °C and 22.3 °C for fluid temperatures of 5 °C, 10 °C and 15 °C respectively. It was interesting to find that although the temperatures at the same location were different for different fluid temperatures, the temperature drop percentages Q were all 27%, which indicated that the influence range of flow channel on heat dissipation was not affected by the fluid temperature Tf.

Influences of Fluid Temperature and Flow Velocity
In order to investigate the influence of fluid temperature on the thermal behavior of the active heat dissipation design, theoretical predictions and FEA were carried out with T f of 5/10/15 • C. The initial and ambient temperature were 25 • C and the fluid velocity was 100 mm/s. Figure 4a shows the temperature along Path 2. Again, there was good agreement between the theoretical prediction and FEA and the maximum temperature difference of 0.56 • C. As expected from FEA, the temperature of the PDMS decreased slightly along the flowing direction, since the heat was taken away by the flowing fluid. This slight temperature change validated the theoretical treatment of the influence of flowing fluid as a forced convection with a constant convective coefficient, although the local convective coefficient decreased along the flowing direction. Figure 4b compares the temperature along Path 3 according to the theoretical prediction and FEA with various fluid temperatures of 5/10/15 • C. It can be observed that the theoretical prediction agreed well with FEA. As expected, lower fluid temperature yielded better heat dissipation performance. At the location of r = 10 mm, the temperatures of PDMS were 19.6 • C, 21 • C and 22.3 • C for fluid temperatures of 5 • C, 10 • C and 15 • C respectively. It was interesting to find that although the temperatures at the same location were different for different fluid temperatures, the temperature drop percentages Q were all 27%, which indicated that the influence range of flow channel on heat dissipation was not affected by the fluid temperature T f .

Influences of Fluid Temperature and Flow Velocity
In order to investigate the influence of fluid temperature on the thermal beha the active heat dissipation design, theoretical predictions and FEA were carried o Tf of 5/10/15 °C. The initial and ambient temperature were 25 °C and the fluid veloc 100 mm/s. Figure 4a shows the temperature along Path 2. Again, there was good ment between the theoretical prediction and FEA and the maximum temperature ence of 0.56 °C. As expected from FEA, the temperature of the PDMS decreased s along the flowing direction, since the heat was taken away by the flowing fluid. Thi temperature change validated the theoretical treatment of the influence of flowin as a forced convection with a constant convective coefficient, although the local con coefficient decreased along the flowing direction. Figure 4b compares the temp along Path 3 according to the theoretical prediction and FEA with various fluid te tures of 5/10/15 °C. It can be observed that the theoretical prediction agreed well wit As expected, lower fluid temperature yielded better heat dissipation performance location of r = 10 mm, the temperatures of PDMS were 19.6 °C, 21 °C and 22.3 °C fo temperatures of 5 °C, 10 °C and 15 °C respectively. It was interesting to find that al the temperatures at the same location were different for different fluid temperatu temperature drop percentages Q were all 27%, which indicated that the influence of flow channel on heat dissipation was not affected by the fluid temperature Tf.   Figure 5a shows the convective coefficient on the flow channel wall as a function of flow velocity from the Sieder-Tate formula. As shown, the convective coefficient increased with increasing fluid velocity. When the flow velocity increased from 20 mm/s to 800 mm/s, the coefficient h f increased from 1172 W/(K·m 2 ) to 4010 W/(K·m 2 ). In addition, the convective coefficient h f increased rapidly when flow velocity v was small, while the changing of coefficient h f became slow when Re exceeded 640. These results indicated that the flow velocity may affect the heat dissipation performance.
the flow velocity may affect the heat dissipation performance.
To demonstrate the influence of flow velocity on the thermal behavior of th heat dissipation design, we show the temperature at the feature point, (Rc, 0) as a of flow velocity v in Figure 5b, where the fluid temperature Tf was 5/10/15 °C. It i that the FEA results were all higher than the theoretical predictions, since the actua local temperatures of the fluid at the fluid-solid interface were accounted for in F responding to a lower value of Tl at the inlet in the theoretical model. When Tf w there was a more distinct difference between the FEA and theoretical prediction. T perature differences were 0.4 °C for the Re of 120 and 0.14 °C for 780. This temp difference decreased with Re because the added heat in the fluid, when flowing f inlet to the feature point (Rc, 0), increased as the flow velocity v slowed, result higher fluid temperature relative to the Tf. Considering that the temperature dif between FEA and theoretical predictions were all lower than 0.4 °C, it can be still ered that the theoretical predictions and FEA were in good agreement.

Thermal Analysis of Active Heat Dissipation Design with Multiple Flow Channels
The previous results and discussion were for the active heat dissipation des a single flow channel. In order to further enhance the heat dissipation performan tiple channels are usually adopted in practical applications. Figure 6a shows a cr tion of the active heat dissipation design with multiple flow channels, with l repr the spacing of neighboring channels. As an analysis of the modelling of multiple c was not possible, only FEA was carried out here to investigate the thermal behav to the symmetric conditions, a single unit cell was developed with a flow channel 1 mm) located in the middle of a square PDMS (in-plane dimension: l  l). The l the flow channel was 60 mm. The bottom and upper sides were set as natural con boundaries with a coefficient of 20 W/(K·m 2 ), and the remaining sides were set as iabatic boundaries. The initial and ambient temperatures were both set at 25 °C. T temperature was 5 °C. To demonstrate the influence of flow velocity on the thermal behavior of the active heat dissipation design, we show the temperature at the feature point, (R c , 0) as a function of flow velocity v in Figure 5b, where the fluid temperature T f was 5/10/15 • C. It is shown that the FEA results were all higher than the theoretical predictions, since the actual higher local temperatures of the fluid at the fluid-solid interface were accounted for in FEA corresponding to a lower value of T l at the inlet in the theoretical model. When T f was 5 • C, there was a more distinct difference between the FEA and theoretical prediction. The temperature differences were 0.4 • C for the Re of 120 and 0.14 • C for 780. This temperature difference decreased with Re because the added heat in the fluid, when flowing from the inlet to the feature point (R c , 0), increased as the flow velocity v slowed, resulting in a higher fluid temperature relative to the T f . Considering that the temperature differences between FEA and theoretical predictions were all lower than 0.4 • C, it can be still considered that the theoretical predictions and FEA were in good agreement.

Thermal Analysis of Active Heat Dissipation Design with Multiple Flow Channels
The previous results and discussion were for the active heat dissipation design with a single flow channel. In order to further enhance the heat dissipation performance, multiple channels are usually adopted in practical applications. Figure 6a shows a cross-section of the active heat dissipation design with multiple flow channels, with l representing the spacing of neighboring channels. As an analysis of the modelling of multiple channels was not possible, only FEA was carried out here to investigate the thermal behavior. Due to the symmetric conditions, a single unit cell was developed with a flow channel (radius: 1 mm) located in the middle of a square PDMS (in-plane dimension: l × l). The length of the flow channel was 60 mm. The bottom and upper sides were set as natural convection boundaries with a coefficient of 20 W/(K·m 2 ), and the remaining sides were set as the adiabatic boundaries. The initial and ambient temperatures were both set at 25 • C. The fluid temperature was 5 • C. Micromachines 2021, 12, x FOR PEER REVIEW 9 of  Figure 6b shows the temperature distribution on the cross-section of the unit cell w a spacing l of 30 mm. s is defined as the distance between the point on the diagonal a the center. The maximum temperature, which occurred at the four corners of the unit ce was 11.2 °C, while the minimum temperature, which occurred at the fluid-solid interfa i.e., where the related temperature drop Q ranged from 69% to 96.5%, was 5.7 °C. Figu 6c shows the temperature distribution along the diagonal with spacing l in the range of to 90 mm. The maximum temperature increased with the increase of spacing, since larger spacing yields a lower heat dissipation capacity, resulting in a higher temperatu When the channel spacing l was 90 mm, the PDMS temperature was 19.6 °C, and Q w only 27% at the position of 15 times the channel radius away from the center compared the same temperature at the position six times the channel radius for the case of a sing flow channel. Thus, the heat dissipation influence of multiple channels is much better th that of a single channel. Figure 6d shows the maximum and minimum temperatures the cross-section with respect to the channel spacing l. The channel spacing had a sign cant influence on the maximum temperature but a negligible influence on the minimu temperature. When the spacing increased from 10 mm to 90 mm, the maximum tempe ture increased by 350%, i.e., from 5.2 °C to 23.4 °C, and Q decreased from 99% to 8%, wh the minimum temperature only increased by 18%, i.e., from 5.0 °C to 5.9 °C, and Q w almost constant. Moreover, when the spacing was larger than the influence range o single channel, its effect became negligible. However, under the same flow rate, t smaller spacing meant that the number of channels increased, which made it necessary improve the pump power to overcome the pressure drop. Therefore, in practical applic tions, the contradiction between the number of channels and cooling performance nee to be adjusted so that cooling capacity and energy saving are maximized.

Conclusions
In summary, a theoretical model based on the method of separation of variables w developed to investigate the thermal behavior of an active heat dissipation design with embedded flow channel. The fluid-solid coupling at the interface was simplified, i.e.,  Figure 6b shows the temperature distribution on the cross-section of the unit cell with a spacing l of 30 mm. s is defined as the distance between the point on the diagonal and the center. The maximum temperature, which occurred at the four corners of the unit cell, was 11.2 • C, while the minimum temperature, which occurred at the fluid-solid interface, i.e., where the related temperature drop Q ranged from 69% to 96.5%, was 5.7 • C. Figure 6c shows the temperature distribution along the diagonal with spacing l in the range of 10 to 90 mm. The maximum temperature increased with the increase of spacing, since a larger spacing yields a lower heat dissipation capacity, resulting in a higher temperature. When the channel spacing l was 90 mm, the PDMS temperature was 19.6 • C, and Q was only 27% at the position of 15 times the channel radius away from the center compared to the same temperature at the position six times the channel radius for the case of a single flow channel. Thus, the heat dissipation influence of multiple channels is much better than that of a single channel. Figure 6d shows the maximum and minimum temperatures on the cross-section with respect to the channel spacing l. The channel spacing had a significant influence on the maximum temperature but a negligible influence on the minimum temperature. When the spacing increased from 10 mm to 90 mm, the maximum temperature increased by 350%, i.e., from 5.2 • C to 23.4 • C, and Q decreased from 99% to 8%, while the minimum temperature only increased by 18%, i.e., from 5.0 • C to 5.9 • C, and Q was almost constant. Moreover, when the spacing was larger than the influence range of a single channel, its effect became negligible. However, under the same flow rate, the smaller spacing meant that the number of channels increased, which made it necessary to improve the pump power to overcome the pressure drop. Therefore, in practical applications, the contradiction between the number of channels and cooling performance needs to be adjusted so that cooling capacity and energy saving are maximized.

Conclusions
In summary, a theoretical model based on the method of separation of variables was developed to investigate the thermal behavior of an active heat dissipation design with an embedded flow channel. The fluid-solid coupling at the interface was simplified, i.e., to forced convention. Finite element analysis and experiments were also carried out to validate the analytical model. It was found that the temperature and velocity of the fluid had distinct heat dissipation effects on the system. Specifically, decreasing the former and increasing the latter significantly improved the heat dissipation performance. Furthermore, the active heat dissipation design with multiple channels was also studied by finite element analysis. It was found that denser flow channels have a much better heat dissipation performance.

Conflicts of Interest:
The authors declare no conflict of interest.