*Article* **Visualization Experimental Study on Silicon-Based Ultra-Thin Loop Heat Pipe Using Deionized Water as Working Fluid**

**Wenzhe Song, Yanfeng Xu, Lihong Xue, Huajie Li and Chunsheng Guo 1, \***

> School of Mechanical, Electrical & Information Engineering, Shandong University, Weihai 264209, China; 201800800305@mail.sdu.edu.cn (W.S.); 201916564@mail.sdu.edu.cn (Y.X.); 202017437@mail.sdu.edu.cn (L.X.); 201800800220@mail.sdu.edu.cn (H.L.)

**\*** Correspondence: guo@sdu.edu.cn; Tel.: +86-631-56887762

**Abstract:** As a type of micro flat loop heat pipe, s-UTLHP (silicon-based ultra-thin loop heat pipe) is of great significance in the field of micro-scale heat dissipation. To prove the feasibility of s-UTLHP with high heat flux in a narrow space, it is necessary to study its heat transfer mechanism visually. In this paper, a structural design of s-UTLHP was proposed, and then, to realize the working fluid charging and visual experiment, an experimental system including a holding module, heating module, cooling module, data acquisition module, and vacuum chamber was proposed. Deionized water was selected as a working fluid in the experiment. The overall and micro phenomena of s-UTLHP during startup, as well as the evaporation and condensation phenomena of s-UTLHP during stable operation, were observed and analyzed. Finally, the failure phenomenon of s-UTLHP was analyzed, and several solutions were proposed. The observed phenomena and experimental conclusions can provide references for further related experimental research.

**Keywords:** loop heat pipe; deionized water; two-phase flow; visualization; heat transfer experiment

#### **1. Introduction**

With the development of electronic technology, as the function of the chip becomes more and more powerful, its power consumption is also rises. At the same time, electronic products are developing in the direction of small, light, and thin, which has caused the heat flux of electronic components to rise sharply. Therefore, the evolution of mobile terminals has led to more stringent requirements on the size of electronic cooling components [1]. Taking smartphones as an example, in recent years, the power consumption of system chips of smartphones has increased to 3–5 W, while the thickness of smartphones has been reduced to about 6 mm [2]. When the equipment runs at high power, it produces a great deal of heat. This heat accumulation leads to uneven temperature distribution, and the resulting thermal stress causes thermal deformation of internal electronic devices [3]. There is evidence indicating that the micro flat loop heat pipe can maximally reduce the temperature of the chip and improve the overall temperature uniformity of the chip. Furthermore, its required heat dissipation space is very small. Thus, the micro flat loop heat pipe appears to be an ideal solution for high-intensity heat dissipation of microscale components [4].

The heat pipe is characterized by high thermal conductivity, long service life, convenient maintenance, and compact and flexible structure. Recently years, research and experiments on heat pipes have attracted more attention [5]. The focus of prior studies mainly differ from the following four aspects: microchannel design, working fluid selection, heat pipe performance, and visualization. In terms of microchannel research, Lim J. [6] proposed a channel layout for flat plate micro heat pipe under local heating conditions that can effectively overcome the limitations of local heating conditions. For working fluid research, Kim J. [7] used ethanol, FC-72, HFE-7000, R-245fa, and R-134a as experimental working fluids to study the selection criteria of working fluids in micro heat pipes. In

**Citation:** Song, W.; Xu, Y.; Xue, L.; Li, H.; Guo, C. Visualization Experimental Study on Silicon-Based Ultra-Thin Loop Heat Pipe Using Deionized Water as Working Fluid. *Micromachines* **2021**, *12*, 1080. https:// doi.org/10.3390/mi12091080

Academic Editors: Junfeng Zhang and Ruijin Wang

Received: 19 July 2021 Accepted: 2 September 2021 Published: 7 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

addition, Narayanasamy M. [8] mixed acetone, deionized water, and tetrahydrofuran fluids with graphene oxide nanoparticles by ultrasonication and prepared nanofluids as working fluids for micro heat pipes for lithium-ion batteries. In terms of performance research, Zhao Y. N. [9] deduced the heat transfer performance of a micro heat pipe array, the influence of inclination angle on heat transfer characteristics, and the pressure limit of the heat pipe structure by using ammonia as a working fluid. Chen G. [10] researched the effects of heat load, cooling water temperature, inclination angle, and other factors on the thermal response performance, temperature distribution, thermal resistance, and other heat transfer characteristics of ultra-thin flat heat pipes. Wang G. [11] concentrated on the effects of working fluid type, charging rate, and inclination angle on the performance of flat plate micro heat pipes through experiments. In terms of visualization research, Kamijima C. [12] studied the relationship between the thermal characteristics and the internal flow characteristics of the micro heat pipe with FC-72 as the working fluid by measuring the effective thermal conductivity, visualizing the flow pattern, and simulating heat transfer according to the flow pattern. The study by Kim Y. B. [13,14] investigated the thermal flow and rapid thermal oscillatory flow in an asymmetric micro pulsation heat exchanger utilizing FC-72 and ethanol as working fluids. It also clarified the optimal charging rate of the two working fluids, evaluated their thermal resistances, identified their key mechanisms of circulating motion, and observed two different flow modes (oscillatory eruption mode and circulating mode). Our research group performed an in-depth study of the characteristics of microchannel flow using numerical simulation [15]. On the basis of our previous research and the four types of research listed above, we conducted in-depth research focusing on visualization and expanded the literature on the design of microchannels and the performance of heat pipes.

Microchannel systems can significantly improve energy transfer efficiency and reduce the overall size of chips. In various technical fields, microchannel technology is one of the most promising areas for the development of equipment [16]. As a kind of phase change heat transfer, gas–liquid two-phase flow heat transfer has very high efficiency [17]. The microchannel is etched on a silicon substrate, and the heat transfer is carried out by gas–liquid two-phase flow, which can realize a flat loop heat pipe with small size and high efficiency. Compared with other working fluids, deionized water has stronger polarity and stronger affinity with silicon chips, which are suitable for silicon-based heat pipes. The phase transition process, flow process, and vapor–liquid distribution of the working fluid can be observed by visualization experiment [18]. Therefore, in this study, a micro ultra-thin loop heat pipe was designed with silicon as the mainboard material, deionized water as the working fluid, and capillary force provided by the microchannel array. The phase transition behavior and vapor–liquid interface of the heat pipe during startup and stable operation were visualized through experiments. The results have an important reference value for the theoretical analysis of the heat transfer mechanism of ultra-thin loop heat pipes.

#### **2. Structure and Experimental Setup of s-UTLHP**

#### *2.1. Structure of s-UTLHP*

The main structure of the s-UTLPH designed in this study, as shown in Figure 1a,b, included a silicon-based motherboard and a Pyrex 7740 glass upper cover plate. The high boron glass upper cover plate was packaged with the silicon-based motherboard by bonding, which allowed us to observe the flow pattern's evolution process in the microchannel visually and intuitively and to understand its physical mechanism. The silicon-based motherboard comprised an evaporation chamber, a condensation chamber, a vapor channel, and a liquid channel connecting the evaporation chamber and the condensation chamber. The overall s-UTLPH size was 40 mm × 20 mm × 1.45 mm. The capillary force of s-UTLPH was provided by parallel microchannels. Both the depth of the microchannel and the whole etching depth were 180 µm. The width of the microchannel was 30 µm, and the aspect ratio was 6:1. In addition, a rectangular vapor overflow cavity was etched on the corresponding

glass cover plate above the parallel microchannel. The height of the vapor overflow cavity was 240 µm, and its structure is displayed in Figure 2.

The evaporation chamber converted the working fluid into vapor by heating in the evaporation chamber, and the capillary core inside provided the capillary driving force to ensure the unidirectional circulation of the working fluid. The main function of the condensing chamber was to condense the vapor working fluid back to the liquid working fluid to release heat. The compensation chamber can store the liquid working fluid of cooling backflow, which was convenient to supplement the liquid working fluid in time and prevent the "dry burning" of the cutoff flow. Eight vapor channels and four liquid channels were located between the evaporation chamber and the condensation chamber. The vapor working fluid produced by phase change flowed to the condensation chamber through the vapor channel. When the condensation chamber became cold, the working fluid turned into liquid and then flowed back to the evaporation chamber through the liquid channel. A hollow insulating slot was also arranged on the s-UTLHP to reduce heat leakage. The size of the s-UTLHP devised in this study is presented in Table 1.

**Figure 1.** s-UTLHP prototype designed in this paper. (**a**) SolidWorks assembly sketch of s-UTLHP; (**b**) Physical drawing of s-UTLHP.

**Figure 2.** Structure diagram of the evaporator.

**Table 1.** s-UTLHP dimension parameters.


#### *2.2. Experimental Device*

To meet the requirements of the working fluid charging, the visualized heat transfer experiment, and the s-UTLHP measurement, it was necessary to arrange the experimental system around the s-UTLHP with limited space. As shown in Figure 3, the visualized heat transfer experimental system involved an s-UTLHP holding module, heating module, cooling module, data acquisition module, vacuum chamber, etc. The data acquisition module consisted of a temperature data acquisition module and a visualization module.

In this study, the silicon-based bottom plate of the s-UTLHP was set with two outlets for charging and pumping. First, a mechanical pump was used for primary air extraction, and an Edward molecular pump was utilized for secondary air extraction. When the pressure was below 1.5 <sup>×</sup> <sup>10</sup>−<sup>3</sup> Pa, the valve at the air extraction port was closed, and then the valve on the other side was opened to fill the working fluid. Deionized water was applied as the working fluid in this paper. The physical parameters of this working fluid are presented in Table 2. After vacuumizing the inner part of the s-UTLHP and its connecting pipelines, the charging operation of the s-UTLHP was carried out. The experimental results demonstrated that high vacuum degree and good charging effect could be obtained by the vacuum charging method, which provides a guarantee for the rapid startup and stable operation of s-UTLHP.

The s-UTLHP, after charging and sealing, was accurately heated by a ceramic heating plate combined with a heat transfer copper block. An ITECH programmable power supply was used of the model IT6121B. Furthermore, a copper cooling block was attached to the condensing end. Constant-temperature cooling water was continuously supplied by a constant-temperature circulating water cooler. The temperature of the condensing end was adjusted by controlling the temperature of the circulating water. In the experiment, the cooling water temperature was set to 5 ◦C. During the experiment, the Omega T thermocouple Was employed to detect the temperature of each key node, and Fluke2638A was used to record and store the temperature data. The arrangement of the thermocouple is displayed in Figure 4. The temperature measuring point was also set to detect the room temperature.

**Figure 3.** Visualization of the heat transfer experiment system.

**Figure 4.** The layout of the thermocouple for temperature measurement. Teva, in—evaporation chamber inlet, Tliq, out—liquid line outlet, Tliq,4,5,6—liquid pipeline, Tc,out—condensation chamber outlet, Tgas—gas channel, Tv—evaporation chamber, Tc,in—condensation chamber inlet.

**Table 2.** Parameters of each liquid working fluid at standard atmospheric pressure and room temperature.


#### *2.3. Error Analysis and Data Processing*

#### 2.3.1. Error Analysis

In this experiment, contact and non-contact temperature measurement methods were adopted through a combination of thermocouple and infrared temperature measurement. In contact measurement, thermocouples were pasted on the s-UTLHP silicon substrate, and thermally conductive silicone grease was coated between the heating block and the silicon substrate to reduce the contact thermal resistance. Meanwhile, the PTFE fixture was designed to reduce the heat loss of the insulation section. The parameters or data to be measured included charge statistics, surface temperature measurement, and input heat calculation. Hence, the main experimental error sources were the temperature measurement point, charge mode, test error, and DC power supply instrument.

(1) Charge uncertainty

A Hamilton 250 µL injector was used to achieve quantitative charging. The minimum calibration interval was 2.5 µL, and the error was ±1.25 µL.

(2) Temperature uncertainty

To ensure the repeatability and effectiveness of the experimental phenomenon, the method of repeated testing was utilized for the experiment under the same conditions. The directly measured parameter value (temperature) in the experiment was calculated according to the following formula:

$$M(T\_{\bar{i}}) = \sqrt{\frac{1}{I(f-1)}\sum\_{j=1}^{J}\sum\_{i=1}^{I}(T\_{\bar{i}j}-\overline{T\_{\bar{i}}})} $$

In the above formula, *J* is the serial number of the repeated experiment and *i* is the serial number of the temperature point tested. In the test results, the maximum deviation of temperature measured by the thermocouple was 1.56 ◦C, and the corresponding maximum uncertainty was 1.93%.

(3) Uncertainty of DC power supply

The uncertainty caused by the instrument is given by class B uncertainty *U* (*E*) = *<sup>A</sup> k* , where A is the nominal error in the instruction manual of the instrument and *k* = √ 3. In the experiment of s-UTLHP heat transfer characteristics, the input heat load was changed by adjusting the voltage of the DC regulated power supply. The nominal error was <sup>±</sup>1%, and the standard uncertainty of voltage was *<sup>U</sup>*(*V*) <sup>=</sup> 1%<sup>√</sup> 3 = 0.577%. According to the principle of error transfer, the uncertainty of input heat load Qin was *U*(*Qin*) = √ (*I*·*u*(*V*))2+(*V*·*u*(*I*))<sup>2</sup> *<sup>V</sup>*·*<sup>I</sup>* <sup>=</sup> 2.341%.

#### 2.3.2. Data Processing

The s-UTLHP charging rate is defined as the ratio of the volume of the working fluid to the total volume of the system. The specific value is computed by the following formula:

$$
\Phi = \frac{V\_{wf}}{V\_{\text{sys}}}
$$

where Φ is the charging rate, *Vw f* is the volume filled with liquid working fluid, and *Vsys* is the volume of the s-UTLHP system.

#### **3. Analysis and Discussion of Experimental Results**

#### *3.1. Startup of the s-UTLHP*

#### 3.1.1. Overall Phenomena at Startup

For the operation of traditional heat pipes, a suitable working fluid should have high latent heat of vaporization and surface tension, low viscosity, and good wetting performance. Additionally, a working fluid with different critical points should be selected to conform to the temperature range [19]. In this study, deionized water was selected as a working fluid.

The startup process of the s-UTLHP filled with deionized water working fluid is presented in Figure 5. The voltage of the regulated power supply was set at 3.4 V, and the current change was recorded after stable operation. The heating power was 5.54 W, the effective heat transfer area was 4 mm <sup>×</sup> 3 mm, and the heat flux was about 46 W/cm<sup>2</sup> .

**Figure 5.** Startup process of the s-UTLHP with two different charging working fluids: (**A**) bubble generation stage at microchannel; (**B**) beginning stage of vapor outflow; (**C**) stable stage after vapor propulsion.

Small bubbles appeared when t = 5 s (Figure 5A). During the experiment, the liquid working fluid in the s-UTLHP evaporation chamber absorbed the external heat through

the chip substrate and increased in temperature. After reaching the saturation temperature, nucleate boiling occurred in the parallel microchannel. The vapor bubbles in the vapor overflow cavity above the microchannel expanded after their appearance to occupy the whole vapor overflow cavity within 5 s (Figure 5B). With a continuous heating process or with an increase in heating power, the vapor bubbles above the vapor overflow cavity decreased and were replaced by pure working fluid vapor. The reason was that the vapor bubbles on the upper cover plate of the evaporation chamber were mixed with the surrounding droplets (liquid film), which was the final result of the dynamic process of evaporation and condensation in the evaporation chamber. Higher heat input made the evaporation effect far greater than the condensation effect, so that purer vapor appeared in the evaporation chamber. In the meantime, the working fluid vapor began to flow from the vapor overflow chamber to the vapor pipe, and the vapor continuously pushed the liquid in the vapor pipe to flow to the condensation chamber (Figure 5C). From Figure 5B,C, it can be seen that the continuous generation of vapor bubbles in the evaporation chamber was weakened, the nucleate boiling was replaced by thin-film evaporation, and the liquid pipe and compensation chamber were filled with liquid and continuously replenished the evaporation chamber. At this time, the s-UTLHP started up and ran stably. In the initial stage of heating and evaporation, it took only 14 s for the deionized water to advance to the vapor channel from the generation of bubbles to the evaporation drying of the vapor overflow chamber.

At the primary stage of startup, the area above the parallel microchannel array of the evaporation chamber evaporated first. During the phase transition of the deionized water working fluid, numerous bubbles were generated, and then a large number of fine droplets were formed, as displayed in Figure 5B,C. Moreover, nucleate boiling emerged in the parallel microchannel at the initial stage of operation.

#### 3.1.2. Micro Phenomena during Startup

To have a more comprehensive understanding of the internal startup and evaporation boiling phenomena for the s-UTLHP, the microscopic phenomena of the s-UTLHP under variable power were observed, and the temperature data during the experiment were recorded. Deionized water has the advantages of being non-toxic and easy to obtain and of having a high latent heat of vaporization, which is commonly used as a working fluid for loop heat pipes [20]. The startup and variable power loading experiments for the s-UTLHP were carried out with deionized water as the working fluid, and the charging rate was 90%. In the experiment, we increased the input power at equal intervals. After each increase of the input power, different heat transfer phenomena occurred. According to these phenomena, we divided the stage into different regions. In addition, each time we heating at the regular rate is a variable power test after reaching steady state at a specific power. The sign of reaching steady state is that the temperature changes very little (gradually flattening from the temperature curve, that is, the temperature is stable or fluctuates within a certain range with the passage of time at a certain power). In this process, the relationship between the overall temperature and power within the experimental time is shown in Figure 6. Ten variable power loading experiments were carried out on the s-UTLHP, and the resulting curve was divided into 11 regions according to the variable power loading time. It can be seen from the figure that the s-UTLHP could respond quickly to different power and heat loads. For the evaporation chamber outlet temperature curve (the highest temperature curve in the figure), after changing the heating power, the temperature first rose rapidly, then gradually decreased and stabilized, and then fluctuated within a certain range. This was because the fluctuation of vapor flow resulted in the fluctuation of heat carried away, which gave rise to the fluctuation of temperature. The temperature of the liquid phase pipeline gradually increased along the fluid return path, because the position near the end of the return pipeline was close to the heating source, which could conduct part of the heat.

According to the above two startup judgment marks for the s-UTLHP, the startup power of the s-UTLHP with deionized water as the working fluid was 3.5517 W, which corresponds to phase E in the curve. The temperature difference between the inlet and outlet of the condensation chamber was 3.18 ◦C. The whole change process of the bubbles in the evaporation chamber during startup (shown in Figure 7) occurred in the area without the shallow cavity and microchannel in the evaporator. The change of bubble volume in stages A–D was not obvious and was accompanied by expansion–contraction oscillation with smaller volume change amplitude, but the magnification was different across the four stages. There were also both large and small droplets, which were constantly generated, merged, grown, and disappeared. The disappearance here was integration either with the liquid around the bubble or with the liquid below. At the E stage, the bubble volume increased rapidly and almost covered the whole area without the shallow cavity and microchannel in the evaporator.

**Figure 6.** Relationship among the temperature and power of s-UTLHP and the time. Teva,in is the inlet temperature of the evaporation chamber; Tliq,out is the outlet temperature of the liquid pipeline; Tliq,4,5,6 is the temperature of the liquid pipeline, evenly arranged along the liquid return direction; Tc,out is the outlet temperature of condensation chamber; Tgas is the temperature of the vapor phase pipeline; Tamb is the ambient temperature; T<sup>v</sup> is the temperature of the evaporation chamber; and Tc,in is the inlet temperature of the condensation chamber.

**Figure 7.** Startup process of the s-UTLHP.

#### *3.2. Stable Operation of the s-UTLHP*

#### 3.2.1. Evaporation Phenomena in Stable Operation

The steady-state performance of the heat pipe was reflected in the variable power operation experiment to a certain extent [21], and the micro phenomena of vapor–liquid distribution, evaporation, and condensation were also captured. In the variable power operation of s-UTLHP, it was found that, as shown in Figure 6, when the power increased, the temperature responded in time, showing a ladder shape on the whole. The temperature curves of the single stages showed a trend of first increasing and then becoming basically stable.

As presented in Figure 8, when deionized water was applied as a working fluid, the s-UTLHP operated stably from stage F. In the subsequent stages, the inlet temperature curve of the vapor phase pipeline was generally stepped, and for each stage, it first increased, then decreased, and then generally stabilized. This was because, for the ceramic heating plate under a given voltage, over time, the current gradually decreased and then stabilized. The microscopic image of stage F at the junction of the evaporation chamber and vapor phase pipeline is presented in Figure 8a. Many droplets can be seen; these droplets constantly repeated the process of generation, growth, disappearance, and washing away. As the power increased, the process changed faster and faster. Figure 8b shows the phenomena in the evaporation chamber. Figure 8c,d show the bubbles at the inlet of the vapor phase pipeline during the stable operation time. When the deionized water working fluid shrank and broke at the inlet of the vapor phase pipeline, slug-like flow appeared, and adjacent bubbles did not easily fuse.

**Figure 8.** Microscopic images of the s-UTLHP at different stages of steady-state operation: (**a**) Image of the connection between the evaporator and the vapor phase line; (**b**) Image inside the evaporator; (**c**,**d**) The image of the inlet of the vapor phase pipeline during stable operation.

#### 3.2.2. Condensation Phenomena in Stable Operation

Like an ordinary loop heat pipe, the s-UTLHP transfers heat through the evaporation– condensation phase transition of working fluid. Figure 9 displays the bubble condensation process in the condensation chamber during stable operation. The positions of Figure 9b–d are illustrated in Figure 9a. These images were obtained under the same magnification and field of view to observe and compare the change of bubble volume. Taking the NCG bubble in focus clearly within the field of view as the initial timing, the 0 s bubble phenomenon is shown in Figure 9b. The condensation process occurred when the heat flux was 37.5575 W/cm<sup>2</sup> .

**Figure 9.** Bubble condensation process in the condensation chamber: (**a**) illustrates where (**b**–**e**) occur in the condenser; (**b**–**e**) are the images at 0 s, 21 s, 45 s and 105 s respectively.

There were many obvious droplets in the deionized water working fluid bubbles. At first, the droplets were small, round, separated from each other, and independent in shape, as shown in Figure 9b. When the edges of the two bubbles grew together as the condensation continued, the condensable components in the bubble gradually condensed. As a result, the droplets were produced or merged with the existing droplets, leading to an increase in the droplet density in the bubble. As shown in Figure 9c,d, the edges of the two bubbles tended to gradually separate. The reason is that the volume of the bubble decreased and the volume of the droplet inside the bubble became larger than depicted in Figure 9c. The two bubbles were completely separated at the time of Figure 9e. The bubble edge of the object was visible, the bubble volume was significantly smaller, and the internal droplet volume became larger.

It was observed that the condensation phenomenon of this s-UTLHP could occur in the gas-phase pipeline when the heat flux was moderate. Bubbles of deionized water working fluid in the gas-phase pipeline during stable operation are shown in Figure 10a–c. The deionized water was ejected into independent bubbles in the gas-phase pipeline, and it was difficult for those bubbles to fuse, even if they shared part of a bubble boundary, as shown in Figure 10a. This was related to the higher surface tension and lower viscosity of deionized water. The volumes of the four bubbles displayed in Figure 10a were different. The lower one had experienced a long time of condensation in the vapor line, while the upper one had just left the column. In the condensation process, slug flow occurred; that is, the vapor column broke at a certain distance, forming independently isolated bubbles. As described above for the condensation chamber, the bubbles were not always completely condensed, and the condensation rate was slow. During operation, the separating vapor

column was positioned at the entrance of the rightmost vapor phase pipeline, as shown in Figure 10c, because of the stress at the corner.

**Figure 10.** Condensation process in vapor phase pipeline.

#### *3.3. Failure of s-UTLHP*

When the heat flux continued to increase, the equilibrium between the refrigerant reflux and phase transition overflow was broken, and the capillary core was considered to be burnt out [22]. Compared with capillary cores made of porous media, the design of a large number of high-aspect-ratio microchannel arrays can make a system have higher stability, effectively reduce the flow resistance, improve the permeability of the evaporation chamber, ensure that the liquid can flow back in time under high heat flux and have higher mass flow rate, and delay the occurrence of flow interruption and drying phenomenon [23]. In this design, the compensation chamber compensated the backflow of liquid in time before the system started up and when the heat flow changed rapidly, so as to improve the stability of the continuous operation of the s-UTLHP. Furthermore, the high-aspect-ratio microchannel array was equipped with a wide pure microchannel area (as shown in area A in Figure 11). This area could generate a great capillary pressure difference through the microchannel under the condition of high heat flux and gas–liquid two-phase transition, thus forming a "hydraulic lock" and effectively preventing backflow. In area B, there was also a microchannel array with a high aspect ratio on the silicon substrate, but unlike in area A, there was a vapor overflow cavity above area B. Figure 11, area C shows the outlet of the vapor overflowing from the evaporation chamber to the gas-phase pipeline, and area D shows the inlet of the liquid working fluid reflux. By setting a reasonable proportion, the circulating pressure change of the working fluid in the s-UTLHP could be effectively balanced.

**Figure 11.** s-UTLHP structure diagram.

An image of the s-UTLHP filled with deionized water working fluid under high heat flux is shown in Figure 12. When the image was taken, the heat input was 41.3 W/cm<sup>2</sup> . It can be observed that under relatively high heat flux, the gas phase pipeline appeared foggy because of two-phase flow. This phenomenon occurred because a great number of gaseous working fluids produced by phase transformation carried small droplets when flowing, and the pressure in the pipeline was high. Furthermore, there was liquid film attached to the wall of the gas phase pipeline. The vapor phase working fluid was produced by continuous phase change in the evaporation chamber. The vapor phase working fluid continuously advanced to occupy most of the gas phase pipeline and was in a relatively stable state.

**Figure 12.** Imaging of the s-UTLHP under high heat flux.

The s-UTLHP in this study had a certain operating limit. When the heat flux reached a certain upper limit, a "dry burning" phenomenon appeared in the evaporation chamber. "Dry burning" refers to the state wherein the working fluid at the inlet and outlet of the evaporation chamber begins to evaporate while the evaporation chamber is continuously heated, finally resulting in the vaporization of all the liquid working fluid in the evaporation chamber. This occurs because the liquid replenishment rate is less than the phase change rate of the working fluid, which means that the liquid cannot be replenished, the liquid film on the wall of the silicon-based microchannel is evaporated quickly, and the gas–liquid two-phase interface regresses until the liquid working fluid in the evaporation chamber completely changes phase. At this time, the single gas-phase state does not have the conditions for capillary force generation and cannot be recycled. As presented in Figure 13, the s-UTLHP filled with deionized water was dried at a heat flux of 44.8 W/cm<sup>2</sup> . Another reason for the drying is that the pressure of the non-condensable gas bubble itself increased due to continuous heating and combined with the vapor. The gas phase flowed back to the evaporation chamber and then dried up instantly; after that, the gas–liquid interface at the gas phase pipeline retrogressed. Because the capillary limit of s-UTLHP was not been reached, the heat flux could further improved if the influence of non-condensable gas was excluded.

In summary, by analyzing the operation state of the s-UTLHP filled with deionized water, we concluded that it was suitable for heat dissipation with high heat flux and environments that are sensitive to temperature fluctuations. In addition, the s-UTLHP failed because of insufficient liquid refrigerant supply after reaching the limit heat flux and because of non-condensable gas reflux.

**Figure 13.** Imaging of s-UTLHP during drying.

#### **4. Conclusions**

In this study, an s-UTLHP with a size of 40 mm in length, 20 mm in width, and 1.5 mm in thickness was designed. The width of the parallel microchannel was 30 µm, and the aspect ratio was 6:1. The design could achieve high capillary force and high mass flow. Furthermore, the longitudinal section of the microchannel with a high aspect ratio could fully extend the liquid film and realize heat exchange with heat flux of 37.5575 W/cm<sup>2</sup> .

Deionized water was selected as the working fluid. In the charging experiment, the s-UTLHP was pumped and charged, and the non-condensable gas was eliminated by combining with visualization. After the s-UTLHP was separately charged with deionized water, it could run stably for a long time. When the heat flux exceeded the limit, the "dry burning" phenomenon was observed, which caused the failure of the s-UTLHP. At low heat flux, the evaporation chamber presented the coexistence of bubbles and droplets, that is, the dynamic equilibrium of evaporation and condensation. At high heat flux, there was pure gas in the evaporation chamber, and the evaporation rate was much higher than the condensation rate, so an efficient evaporation replenishment mechanism was carried out.

This study explored a self-circulating heat dissipation scheme suitable for microchip integration, realized phase-change heat transfer on the surface of a silicon substrate, and proved the feasibility of the scheme for high-heat-flux heat dissipation in a narrow space. Moreover, the regular structure of the high-aspect-ratio microchannel array in the core of the evaporation chamber means that these arrays can be mass-produced by a 3D IC manufacturing process, which could help with popularization and application. In addition, this research adopted a visual method to conduct macroscopic and microscopic observational studies of evaporation and condensation, respectively, and to further explore the boiling flow and condensation cycle characteristics of the working fluid in the microchannel. It also provided references for the improvement of fine structure and surface morphology and for the selection of a working fluid, which can enhance the heat transfer capacity and the stability of the cyclic heat transfer mechanism. Therefore, this study has important reference value in the field of micro-scale heat dissipation for electronic chips.

**Author Contributions:** Conceptualization, C.G.; data curation, W.S.; formal analysis, W.S.; funding acquisition, C.G.; investigation, H.L.; methodology, Y.X.; project administration, C.G.; software, L.X.; supervision, C.G.; validation, Y.X.; visualization, W.S. and Y.X.; writing—original draft, W.S.; writing—review and editing, W.S., C.G., and Y.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Focus on Research and Development Plan in Shandong Province (Grant No. 2018GGX104011) and the Natural Science Foundation of Shandong Province (Grant No. ZR2017BEE012).

**Data Availability Statement:** The data that support the findings of this study are available from the authors upon reasonable request. The data are not publicly available due to that they are planned to be used in our following experiments.

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

#### **Nomenclature**


#### **References**


### *Article* **Effects of Microscopic Properties on Macroscopic Thermal Conductivity for Convective Heat Transfer in Porous Materials**

**Mayssaa Jbeili and Junfeng Zhang \***

Bharti School of Engineering and Computer Science, Laurentian University, 935 Ramsey Lake Road, Sudbury, ON P3E 2C6, Canada; mjbeili@laurentian.ca

**\*** Correspondence: jzhang@laurentian.ca; Tel.: +1-705-675-1151 (ext. 2248)

**Abstract:** Porous materials are widely used in many heat transfer applications. Modeling porous materials at the microscopic level can accurately incorporate the detailed structure and substance parameters and thus provides valuable information for the complex heat transfer processes in such media. In this study, we use the generalized periodic boundary condition for pore-scale simulations of thermal flows in porous materials. A two-dimensional porous model consisting of circular solid domains is considered, and comprehensive simulations are performed to study the influences on macroscopic thermal conductivity from several microscopic system parameters, including the porosity, Reynolds number, and periodic unit aspect ratio and the thermal conductance at the solid–fluid interface. Our results show that, even at the same porosity and Reynolds number, the aspect ratio of the periodic unit and the interfacial thermal conductance can significantly affect the macroscopic thermal behaviors of porous materials. Qualitative analysis is also provided to relate the apparent thermal conductivity to the complex flow and temperature distributions in the microscopic porous structure. The method, findings and discussions presented in this paper could be useful for fundamental studies, material development, and engineering applications of porous thermal flow systems.

**Keywords:** heat transfer; porous media; pore-scale modeling; boundary condition; thermal conductivity; porosity; conjugate interface; aspect ratio

#### **1. Introduction**

Fluid flows and the associated heat transfer in porous materials have attracted great interest over the decades for their scientific and practical importance. With the large contact area between the solid matrix and the coolant fluid, heat transfer performance is significantly enhanced in such porous materials. Relevant applications can be found in various industrial areas, including catalytic beds for chemical reactions, water treatment, nanofluids, heat exchangers, heat sinks and thermal energy storage systems [1–8]. Extensive experimental and theoretical investigations have been conducted to study the heat transfer performance in porous materials [9–11]. In addition to the high cost and long time duration, experimental studies are not able to reveal the flow and temperature distributions in the complex microscopic porous geometry, while such information is crucial for understanding the mechanism and relationships among various factors involved. On the other hand, theoretical analysis is limited to simple geometry structures and idealized flow situations, and such conditions can hardly be satisfied in realistic systems. Fortunately, with the advances in computer and software technologies, numerical simulations have been proved to be useful for various complex systems, including flows and heat transfer in porous media [12–16]. In the literature, there exist two typical approaches for modeling the flow and heat transfer in porous media: the continuum and the pore-scale approaches. The first method treats the porous materials as continuum media and solves the macroscopic flow and energy equations with apparent parameters [17,18]. Despite the model's simplicity and computational efficiency, this approach relies on the accuracy of those apparent

**Citation:** Jbeili, M.; Zhang, J. Effects of Microscopic Properties on Macroscopic Thermal Conductivity for Convective Heat Transfer in Porous Materials. *Micromachines* **2021**, *12*, 1369. https://doi.org/10.3390/ mi12111369

Academic Editor: Kwang-Yong Kim

Received: 11 October 2021 Accepted: 5 November 2021 Published: 7 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

parameters and also cannot incorporate the specific microscopic structure and material thermophysical properties in the porous media. On the other hand, the pore-scale method solves the flow and energy equations with the microscopic structure and properties of the matrix material considered explicitly [19,20]. Pore-scale simulations can provide detailed flow and thermal fields at the microscopic level, and apparent properties can be obtained from these microscopic distributions for macroscopic analysis [21–23]. In recent years, the lattice Boltzmann method (LBM) has become especially popular in such simulations, mainly for its convenience in dealing with complex boundary geometry [22–25].

With the microscopic porous structure modeled explicitly, the physical scale of porescale simulations is limited to small sizes due to the large computation demand. For this reason, pore-scale simulations typically work with a small unit and assume that such identical units are repeated in space. To solve the flow and energy equations in a periodic unit, appropriate boundary conditions are necessary on the side surfaces of this computational unit. Such side surfaces are simply for computational convenience and they are not physical surfaces or boundaries; therefore, definite conditions on such surfaces are not available. To create a thermal gradient, previous studies are usually assigned constant temperatures on two opposite unit surfaces and assumed adiabatic (no heat flux) conditions on other side walls [24–27]. This treatment has been examined recently and significant concerns have been raised in terms of its validity and accuracy [28]. Few studies treated the unit surface conditions more reasonably by considering the periodic relations of temperature, however, the temperature or heat flux on the fluid-solid interface was fixed at constant values [29–31]. Clearly this is physically not true: the thermal condition (temperature and heat flux) at the microscopic fluid-solid interface is determined by the particular local flow and thermal situations and it cannot be specified in advance. The interfacial thermal resistance may also exist at the fluid-solid interface in porous materials [32] or the constituent interface in composite materials [33], and this feature has typically been neglected in previous studies. Furthermore, previous studies usually only considered situations with the flow along the temperature gradient direction [22,25,30], whereas in many experimental setups and industrial applications the temperature gradient is mainly in the orthogonal plane to the fluid flow [34–36].

In this paper, we study the convective heat transfer in porous materials using the recent generalized periodic boundary method for thermal flows by Jbeili and Zhang [28]. This boundary method was established with rigorous mathematical justifications and validated by carefully designed numerical tests. Using this method, extensive simulations are conducted based on a two-dimensional (2D) porous model of circular solid particles. Unlike previous studies usually focusing on the porosity and Reynolds number effects, we also investigate the influences of the unit aspect ratio and the interface conductance on the macroscopic thermal performances of porous flows. Our results show that, even with the same porosity and Reynolds number, the unit aspect ratio can greatly affect the apparent thermal conductivity at the microscopic level. It is also interesting to observe that a weaker interfacial conductance can reduce the tortuosity conductivity but increase the dispersion conductivity, and that the dispersion conductivity can respond to the porosity change in a non-monotonic fashion. In addition, qualitative discussions are provided to relate the macroscopic conductivity coefficients to the microscopic flow and thermal fields for a better understanding of the complex nature of porous thermal flows.

#### **2. Governing Equations and Boundary Conditions for Periodic Unit**

For simplicity and clarity, in this study we consider the 2D porous model with circular solid domains as shown in Figure 1a. The fluid flows vertically from the bottom to the top, while the thermal gradient is imposed in the horizontal direction. The governing equations for this thermal flow system include the continuity equation, Equation (1), and the Navier–Stokes equation Equation (2) for the fluid domain Ω0:

$$\nabla \cdot \mathbf{u} = \mathbf{0},\tag{1}$$

$$
\rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p + \mu \nabla^2 \mathbf{u} + \mathbf{f};\tag{2}
$$

and the energy equation

$$\frac{\partial T\_f}{\partial t} + \mathbf{u} \cdot \nabla T\_f = \mathfrak{a}\_f \nabla^2 T\_{f\prime} \quad \text{in } \Omega\_{0\prime}$$

$$\frac{\partial T\_s}{\partial t} = \mathfrak{a}\_s \nabla^2 T\_{s\prime} \quad \text{in } \Omega\_{1\prime} \tag{3}$$

for the fluid (Ω0) and solid (Ω1) domains, respectively [37]. In these equations, *ρ* represents the fluid density, **u** is the fluid velocity, *p* is the pressure, *µ* is the dynamic viscosity of fluid, **f** is the body force on the fluid, *T* is the temperature, and *α* is the thermal diffusivity. The subscripts *<sup>f</sup>* and *<sup>s</sup>* for *T* and *α* denote the fluid and solid domains respectively. On the fluid-solid interface Γ, the no-slip boundary condition is applied for the fluid flow

$$\mathbf{u} = \mathbf{0}, \quad \text{on} \quad \Gamma,\tag{4}$$

to incorporate the possible thermal resistance and therefore the temperature discontinuity at the solid–fluid interface Γ, a general conjugate condition for temperature is given as [38–40]:

$$q = \mathcal{C} \left( T\_f - T\_s \right) = K\_s \frac{\partial T\_s}{\partial n} = K\_f \frac{\partial T\_f}{\partial n}, \quad \text{on } \Gamma. \tag{5}$$

Here, *q* is the heat flux along the local normal direction **n**, which points from Domain Ω<sup>1</sup> toward Domain Ω<sup>0</sup> (Figure 1). *C* is the interface conductance, and *K<sup>f</sup>* and *K<sup>s</sup>* are the thermal conductivities of the fluid and solid substances, respectively.

**Figure 1.** Schematic representations of the 2D porous model used in this study (**a**) and one periodic unit considered in pore-scale simulations (the dashed box in (**a**), and with more details in (**b**). The key parameters involved in the model description are also labeled, and more details can be found in the text.

In the pore-scale approach, only one periodic unit of the porous material is considered in the simulation (Figure 1). To solve the governing equations in this periodic unit, appropriate boundary conditions are required for the side boundaries of the periodic unit. These side boundaries are virtual but not real surfaces, and thus physical constraints on such boundaries are not directly available. For flows in such periodic structures with negligible fluid property changes, it has been well recognized that the flow field would be identical in each periodic unit [21,41–43]:

$$\mathbf{u}(\mathbf{x} + ml, y + nl) = \mathbf{u}(\mathbf{x}, y),\tag{6}$$

where *l* and *h* are the unit length and height (Figure 1), and *m* and *n* are arbitrary integers. Based on previous practices in periodic thermal flow simulations [28,41–45], the following relationship is proposed for the temperature field among periodic units in Figure 1a:

$$T(\mathbf{x} + ml, y + nh) = T(\mathbf{x}, y) + mlT\_{\mathcal{g}}'.\tag{7}$$

Here, *T* ′ *g* is the global thermal gradient in the horizontal direction, and it can be obtained from the temperature difference *T<sup>H</sup>* − *T<sup>L</sup>* and the distance *L* between the locations where these two temperatures *T<sup>L</sup>* and *T<sup>H</sup>* are imposed (see Figure 1a) as *T* ′ *<sup>g</sup>* = (*T<sup>H</sup>* − *TL*)/*L*.

Assume **u**(*x*, *y*) and *T*(*x*, *y*) are the correct flow and temperature solutions in one periodic unit (0 ≤ *x* ≤ *l* and 0 ≤ *y* ≤ *h*), it can be readily shown that the velocity from Equation (6) and temperature from Equation (7) satisfy all the governing equations and boundary conditions given in Equations (1)–(5), and thus they are the correct solutions for the unit of *ml* ≤ *x* ≤ (*m* + 1)*l* and *nh* ≤ *y* ≤ (*n* + 1)*h*. Therefore, in pore-scale simulations, one can simulate the flow and thermal fields in one periodic unit as shown in Figure 1b, and the solutions can be extended to other units according to Equations (6) and (7). These relations can then be used to establish correct boundary conditions for the side surfaces of the periodic unit. For example, the right boundary *x* = *l* for the unit shown in Figure 1b actually is also the left boundary of the next unit on its right. According to Equations (6) and (7) with *m* = 1 and *n* = 0, we have the left-right boundary relations for the periodic unit as:

$$\mathbf{u}(l,y) = \mathbf{u}(0,y), \quad T(l,y) = T(0,y) + lT'\_{\mathcal{S}}.\tag{8}$$

Similarly, the top-bottom boundary relations for flow and temperature are expressed as:

$$\mathbf{u}(\mathbf{x},h) = \mathbf{u}(\mathbf{x},0), \quad T(\mathbf{x},h) = T(\mathbf{x},0). \tag{9}$$

These relations are similar to those used by Kuwahara et al. [21]; however, no mathematical justifications and numerical validations were provided there. In addition to these boundary relations, due to the linearity and homogeneity of the governing equations and boundary conditions, extra anchoring conditions are necessary to ensure the numerical convergence [39,41,42]. For our current system in Figure 1a, the following conditions are adopted in our next simulations:

$$\frac{1}{l}\int\_{0}^{l}u(\mathbf{x},0)d\mathbf{x} = \mathcal{U}\_{0}, \quad \frac{1}{h}\int\_{0}^{h}T(0,y)dy = T\_{0};\tag{10}$$

where *U*<sup>0</sup> is the mean velocity through the porous medium and *T*<sup>0</sup> is the mean temperature at the left unit boundary. The mean velocity *U*<sup>0</sup> is related to the Reynolds number *Re* of the system (to be defined later), while the mean temperature *T*<sup>0</sup> can be set at an arbitrary value and it has no impact on the result analysis and interpretation.

#### **3. Numerical Validation of the Periodic Relations among Units**

To further verify the flow and temperature relations among units given in Equations (6) and (7), we conduct a direct numerical simulation for the system in Figure 1a. This simulation is also helpful to illustrate our concerns with the periodic unit boundary treatments used in previous studies. The diameter of the solid domains is *D* = 50 and they are arranged as a square array with *l* = *h* = 2*D*. The porosity of this model medium is thus *ǫ* = *π*/16. The rectangular simulation domain has a length of *L* = 20*D* in the horizontal direction and a hight of *H* = 10*D* in the vertical direction, and thus the simulation domain consists 10 × 5 = 50 identical periodic units as shown in Figure 1b. Constant temperatures are imposed on the left and right boundaries: *T<sup>L</sup>* = 0 at *x* = 0 and *T<sup>H</sup>* = 1 at *x* = *L*, resulting a global thermal gradient of *T* ′ *<sup>g</sup>* = (*T<sup>H</sup>* <sup>−</sup> *<sup>T</sup>L*)/*<sup>L</sup>* = 10−<sup>3</sup> . The classical periodic condition is imposed on the top and bottom boundaries for both flow velocity and temperature, and on the left and right boundaries for the fluid flow:

**u**(0, *y*) = **u**(*L*, *y*), **u**(*x*, 0) = **u**(*x*, *H*), and *T*(*x*, 0) = *T*(*x*, *H*). Please note that the temperature relation Equation (7) is not involved in this validation simulation. A body force **f** in the upward direction is utilized to generate the flow, and its magnitude is adjusted to obtain the Reynolds number *Re* = *ρU*<sup>0</sup> √ *lh*/*µ* = 50. Here we use the geometric mean value of the periodic unit length *l* and height *h* for the Reynolds number *Re*, and this choice is made for our convenience in examining the effect of the aspect ratio *β* = *h*/*l* in Section 4.1. The Prandtl number of the fluid is 0.7. The thermal conductivities are set as *K<sup>f</sup>* = 0.2 for fluid and *K<sup>s</sup>* = 10*K<sup>f</sup>* = 2 for solid. A relatively large interface conductance value *C* = 10 is adopted to minimize the temperature discontinuity across the interface. The lattice Boltzmann method (LBM) with the D2Q9 (2D and with nine lattice velocities) lattice structure is used to solve the flow and thermal fields [44,46–48] for all calculations in this paper, and the counter-extrapolation method [49] is adopted for the conjugate thermal condition on domain interfaces. As in general LBM studies [16,22,31,48–50], all quantities provided above and in the later result presentation are non-dimensional based on LBM simulation units (for example, length in the lattice grid resolution, and time in the simulation time step).

The calculated flow and temperature fields from this direct simulation are displayed in Figure 2. The flow pattern in Figure 2a appears identical in each periodic unit. For the temperature in Figure 2b, it can be seen that the temperature increases in general along the horizontal direction, agreeing with the global thermal gradient generated from the two boundary temperatures *T<sup>L</sup>* = 0 at the left and *T<sup>H</sup>* = 1 at the right. In Figure 2b, neither the temperature nor the heat flux are constant along the solid–fluid interface (edges of the circular patches), meaning that the constant temperature or flux assumptions for the interfaces in previous studies [22,30,31] are not valid. The wider gaps between isotherm contours in the solid domain imply smaller temperature gradients there, and this is consistent with the larger solid conductivity used in this simulation. With a relatively larger conductance value for the solid–fluid interface, no apparent discontinuity in temperature can be observed across the interfaces.

**Figure 2.** Flow (**a**) and temperature (**b**) distributions obtained from the direct simulation for the 2D porous model system in Figure 1a.

The apparently linear increase in temperature and the similarity in local isotherm contours in Figure 2b appear rational according to the temperature relation Equation (7). Moreover, *T* ∼ *y* profiles at six horizontal locations are plotted in Figure 3, one in each panel from the left and with the *x* positions labeled on top. These horizontal positions are selected at the same relative locations to the nearby solid columns (0.4*D* to the left of the patch centers). The flat segments on these curves occur in the solid domains, and they are due to the vertical isotherm lines there. The strong variations in these profiles clearly show that in general one should not assign constant temperatures on the periodic unit boundaries, as done in [24–26]. According to Equation (7), these temperature profiles should have the same variation features along the vertical direction, except different offset values. This is apparently true by looking at these individual profiles. For a more quantitative

confirmation, we then shift each profile by its mean temperature *Tm*, and the six shifted profiles are plotted in the last panel on the right in Figure 3. These six shifted profiles overlap with each other perfectly and we only see one single curve there. These identical shifted profiles precisely confirm that the temperature values at same relative positions in individual periodic units are only different in respective mean temperatures and the variation fashion is identical. Furthermore, the mean temperature values calculated along these six profiles are listed in Table 1. According to Equation (7), these mean temperatures can be related to the mean value at *x*/*D* = 0.6 by

$$T\_m(\mathbf{x}\_i) = T\_m(\mathbf{x} = 0.6D) + (\mathbf{x}\_i - 0.6D)T\_{\mathcal{S}'}' \tag{11}$$

with *x<sup>i</sup>* = 0.6*D*, 2.6*D*, 6.6*D*, 10.6*D*, 14.6*D* and 18.6*D* for the six profiles in Figure 3 (from left to right). The predicted values from this equation are also provided in Table 1 for comparison, and the excellent agreement convincingly confirms the validity of the temperature relation Equation (7) for such porous thermal flows.

**Figure 3.** Temperature profiles along the vertical direction at different horizontal positions (labeled on top) from the direct simulation results in Figure 2b. The last panel on the right collects all the temperature profiles in other panels, however shifted by their individual mean temperatures. These six shifted profiles become identical and completely overlap each other, confirming the thermal relation Equation (7) among periodic units along the horizontal direction.

**Table 1.** Comparison of the mean temperature values at different horizontal positions obtained from the simulation and predicted by Equation (11).


#### **4. Simulation Results and Discussion**

The enhanced heat transfer performance in porous flow systems, in principle, benefits from the large solid–fluid contact area and the long, twisted flow path as the fluid passes through the microscopic porous structure. At the macroscopic, practical level, the apparent thermal conductivity is often used to quantify the overall thermal behaviors. The volume averaging analysis [21,51–53] can be applied to obtain the effective conductive tensor from the microscopic flow and thermal fields. Kuwahara et al. [21] further decomposed this conductivity tensor **K** into three parts

$$\mathbf{K} = \mathbf{K}\_{\ell}\mathbf{I} + \mathbf{K}\_{tor} + \mathbf{K}\_{dis\prime} \tag{12}$$

where *Ke* is the stagnant conductivity, which is simply the volume average of the fluid and solid conductivities based on the porosity *ǫ*:

$$\mathbf{K}\_{\ell} = \varepsilon \mathbf{K}\_{f} + (1 - \varepsilon)\mathbf{K}\_{\text{s}}.\tag{13}$$

Please note that the expression for *Ke* in Equation (13) is simply the coefficient for the isotropic part of the conductivity matrix **K** from the volume average analysis [21,51–53]; and that there are no assumptions involved on the microscopic porous structures. **K***tor* and **K***dis* are called the tortuosity and dispersion conductivity tensors, respectively. When the mean flow direction is along one of the coordinate directions, only diagonal elements are nonzero; and obviously the primary concerns are the diagonal elements in the global thermal gradient direction [21,53,54]. For the setup in Figure 1, we simply use *Ktor* and *Kdis* to denote the *xx* components of tensors **K***tor* and **K***dis*; and they can be calculated from the simulated flow and temperature distributions in one periodic unit by [21,53]

$$K\_{tor} = -\frac{K\_{\rm s} - K\_f}{h l T\_{\rm g}'} \int\_{\Gamma} n\_x T\_f d\Gamma\_{\prime} \tag{14}$$

$$K\_{\rm dis} = -\frac{\rho c\_{pf}}{h l T\_{\rm g}'} \int\_{\Omega\_0} (T - \vec{T}) (\mathbf{u} - \vec{\mathbf{u}}\_f)\_{\rm x} d\Omega;\tag{15}$$

where *cp f* is the fluid specific heat, and the volume average temperature *T*¯ and the intrinsic average velocity **u**¯ *<sup>f</sup>* are given by:

$$\hat{T} = \frac{1}{l\hbar} \left( \int\_{\Omega\_0} T\_f d\Omega + \int\_{\Omega\_1} T\_s d\Omega \right), \quad \hat{\mathbf{u}}\_f = \frac{1}{\epsilon l\hbar} \int\_{\Omega\_0} \mathbf{u} d\Omega. \tag{16}$$

The subscript *x* denotes the *x* component of the corresponding vector.

In this section, we apply the boundary conditions in Section 2 to the periodic unit shown in Figure 1b, and examine the effects on the effective macroscopic conductivity coefficients *Ktor* and *Kdis* from several microscopic parameters, including the aspect ratio *β* = *h*/*l*, the Reynolds number *Re* = *ρ* √ *hlU*0/*µ*, the interfacial thermal conductance *C*, and the porosity *<sup>ǫ</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>π</sup>D*2/4*hl*. We keep the unit area *lh* <sup>=</sup> 90, 000 constant in our next simulations, so that the Reynolds number *Re* is directly proportional to the mean flow velocity *U*0, and the porosity *ǫ* depends purely on the solid diameter *D* and not affected by the aspect ratio *β*. The mean temperature at *x* = 0 is set as *T*<sup>0</sup> = −0.05 and the global thermal gradient is *T* ′ *<sup>g</sup>* = 0.1/*l*. Please note that, due to the linearity and homogeneity of the energy equation Equation (3) and the associated thermal boundary conditions, the particular values of *T*<sup>0</sup> and *T* ′ *<sup>g</sup>* used in the simulations do not affect the calculated conductivities *Ktor* and *Kdis* from Equations (14) and (15). To make this point clear, we consider a periodic unit under two situations: Situation (a) with mean inlet temperature *T*0,*<sup>a</sup>* and thermal gradient *T* ′ *g*,*a* ; and Situation (b) with mean inlet temperature *T*0,*<sup>b</sup>* and thermal gradient *T* ′ *g*,*b* . If *Ta*(*x*, *y*) is the solution in the periodic unit with Situation (a), the temperature field for Situation (b) should be *T<sup>b</sup>* (*x*, *y*) = *T*0,*<sup>b</sup>* + *T* ′ *g*,*b T* ′ *g*,*a* [*Ta*(*x*, *y*) − *T*0,*a*]. Clearly, *T<sup>b</sup>* (*x*, *y*) satisfies the energy equation Equation (3) and all temperature conditions in Equations(5), (8), (9), (10), . Substituting this expression of *T<sup>b</sup>* (*x*, *y*) in Equations (14) and (15) and considering R Γ *uxd*Γ = 0 for periodic structures, one can find that the *Ktor* and *Kdis* values remain the same for Situations (a) and (b). Other simulation parameters are kept the same as given in Section 3, unless otherwise mentioned. As for the initial state, we start simulations with a linear

transition from *T*<sup>0</sup> at the inlet to *T*<sup>0</sup> + *T* ′ *g l* at the outlet for the temperature, and zero velocity for the flow. The results are taken when steady or quasi-steady states are established. The computer code for this work has been developed based on our programs used in previous publications [28,38,44,49,55,56], and all important elements involved in our simulations, including the LBM algorithms for flow and heat transfer, the no-slip boundary and conjugate interface treatments, as well as the mesh resolution selection, have been validated and confirmed in these previous studies.

#### *4.1. Effects of the Aspect Ratio β and Reynolds Number Re*

We start with simulations to study the effect of the aspect ratio *β* on the macroscopic thermal conductivity, which has not been well addressed in previous investigations. Here we set the porosity *ǫ* = 0.85, and accordingly the solid domain diameter is *D* = 131.11. Following previous studies [29,57,58], two representative Reynolds number values, *Re* = 50 and 100, are considered in this section. Higher Reynolds numbers are possible for gas flows through porous media [59,60]. Our simulations cover a range of 0.25∼4 for the aspect ratio *β*; further increasing (>4) or decreasing (<0.25) of the *β* value will make the gap between the solid surfaces too small for accuracy and stable computations. The flows are always steady for all tested *β* values at *Re* = 50 and for *β* ≤ 1 at *Re* = 100; however, the flow becomes unsteady for *β* > 1 and *Re* = 100. This transition is reasonable, since for a larger *β*, the gap between the solid surfaces becomes narrower, and the fluid velocity through this gap increases accordingly. This situation is similar to a flow jet coming out from a small opening shooting into a relatively large space. Compared to the condition with the same mean flow velocity (same Reynolds number) but a smaller aspect ratio *β*, the reduction in flow passage width is less significant (a wider gap between solid surfaces) and the fluid space after the gap is relatively limited. Figure 4 shows the streamlines and temperature distributions for two representative aspect ratios, *β* = 0.44 and 1.44 at *Re* = 50 and 100. For steady flows in Figure 4a–c, we see smooth and symmetric streamlines and isotherm contours. On the other hand, when the flow become unsteady in Figure 4d, these lines are distorted and less organized. Both for the steady and unsteady cases, thermal features discussed in Section 3 are noticed as well, including the relatively uniform temperature distributions inside the solid domains (due to the high solid conductivity *Ks*) and the smooth temperature transition across the solid–fluid interface (due to the large interface conductance *C*). Moreover, we see significant temperature variations along the unit boundaries for the unsteady system in Figure 4d. The variation fashions are similar on opposite edges. More specific, the temperature variations along the top and bottom boundaries appear identical, whereas the temperature is low on the left and high on the right. All these observations are consistent with the thermal boundary conditions in Equations (8) and (9).

For unsteady flow situations, the mean flow velocity *U*<sup>0</sup> and Reynolds number *Re* both vary with time, and it is difficult to achieve the time-averaged Reynolds number exactly at 100. For such situations, we accept the simulation results when the *Re* variation is limited in the range of 95 ∼ 105. The unsteady flow certainly affects the thermal field, and accordingly, the tortuosity and dispersion conductivities *Ktor* and *Kdis* calculated from Equations (14) and (15) also vary with time. Figure 5 displays the temporal oscillations of *Re*, *Ktor* and *Kdis* for *Re* = 100 and *β* = 1.44. The oscillation periods for *Ktor* and *Kdis* appear to be identical, as twice that for *Re*. In our next analysis, we use the average values over a variation period for the tortuosity and dispersion conductivities *Ktor* and *Kdis*.

Figure 6 plots the apparent conductivities *Ktor* and *Kdis* changing with the aspect ratio *β* and Reynolds number *Re*. Clearly, even under the same porosity *ǫ* = 0.85, the macroscopic conductivities can be greatly affected by the aspect ratio *β*, and the influence can be enhanced with a higher Reynolds number *Re*. The tortuosity conductivity *Ktor* decreases with *β* in Figure 6a, which can be explained by looking at the *Ktor* definition Equation (14) and the temperature distributions in Figure 4. Equation (14) can be interpreted as a weighted sum of the fluid temperature over the solid–fluid interface Γ, with *nx*, the *x*-component

of the normal vector **n**, as the weighting factor. For the system considered here, the left surface has *n<sup>x</sup>* > 0 and the right surface has *n<sup>x</sup>* < 0. On each semicircular surface, the portion around the horizontal centerline (the surface is approximately aligned in the vertical direction and thus it has a larger |*nx*| value) plays a more determinant role than the portions near the unit edges (the surface is approximately aligned in the horizontal and thus |*nx*| ≈ 0). Therefore, by comparing the temperatures on the parts around the horizontal centerline of the two semicircular surfaces in Figure 4, we can have a qualitative understanding of the *Ktor* dependence on *β*: For a large *β*, the two semicircular surfaces are closer and the temperature difference on the closest portions becomes smaller, and this smaller temperature difference results in a smaller tortuosity conductivity, as observed in Figure 6a. The Reynolds number *Re* appears to be less influential on *Ktor*, except for the highly unsteady case with *β* = 4. The weak influence of *Re* on *Ktor* is consistent with that observed in Ref. [21].

**Figure 4.** Simulation results of the flow streamlines and temperature distributions for *β* = 0.44 (**a**,**b**) and *β* = 1.44 (**c**,**d**) at *Re* = 50 (**a**,**c**) and 100 (**b**,**d**).

Unlike the tortuosity conductivity *Ktor*, the dispersion conductivity *Kdis* increases by orders with *β* in Figure 6b. Similarly, we attempt to interpret this observation by examining the *Kdis* definition Equation (15) and flow and thermal fields in Figure 4. Clearly, Equation (15) represents the disorder level of the flow and thermal distributions in the periodic unit, since the integrand is simply the product of the temperature fluctuation *<sup>T</sup>* <sup>−</sup> *<sup>T</sup>*¯ and the flow fluctuation (**<sup>u</sup>** <sup>−</sup> **<sup>u</sup>**¯ *<sup>f</sup>* )*x*. Therefore, the more disordered the flow and temperature distributions in the periodic unit, the larger *Kdis* value will be obtained. This analysis is consistent with our general understanding of the thermal dispersion process and it explains the increasing trend of *Kdis* with *β* and *Re* in Figure 6b as well. The similarity in flow and thermal fields between *Re* = 50 and 100 for *β* ≤ 1 steady systems suggests that *Kdis* does not change much with *Re* there; however, for *β* > 1, the *Re* = 100 systems become unsteady and the flow and temperature distributions become significantly disordered, resulting in an abrupt increase in *Kdis*. We also notice that only the *x*-component of the

flow fluctuation is involved in Equation (15). For the steady flows at *β* ≤ 1, the nonzero *x* fluid velocity is limited in the small circulation areas near the four unit corners, whereas for large *β* values, the circulation region is large. The high Reynolds number *Re* = 100 is also helpful in increasing the circulation size, and it can even make the *x* velocity component nonzero over almost the entire fluid domain for the unsteady flows with *β* > 1. This is the reason that we see the *Kdis* increase with *β* and *Re* by orders for *β* > 1 in Figure 6b.

**Figure 5.** Variations of *Re* (**a**), *Ktor* (**b**) and *Kdis* (**c**) with time for the unsteady system with *β* = 1.44 and *Re* ≈ 100. The simulation time is normalized based on the mean velocity *U*<sup>0</sup> and solid cylinder diameter *D*.

**Figure 6.** Plots of *Ktor* (**a**) and *Kdis* (**b**) changing with the aspect ratio *β* at two Reynolds numbers *Re* = 50 (black squares) and 100 (blue circles).

#### *4.2. Effect of the Interfacial Conductance C*

Next we shift our attention to the influence of the interfacial thermal conductance *C* on the apparent conductivities. Again this is an important topic [32,33], however, it has not been investigated adequately. In this part, we fix the porosity *ǫ* = 0.85 and the Reynolds number *Re* = 50 for simplicity, and test three interfacial conductance values *C* = 10 (virtually no thermal resistance at the solid–fluid interface as observed in Figures 2b and 4), <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> , and 5 <sup>×</sup> <sup>10</sup>−<sup>5</sup> . Figure 7 displays the flow and thermal distributions at these three conductance values with two aspect ratios *β* = 0.44 and 1.44; and Figure 8 collects the calculated *Ktor* and *Kdis* conductivities in our simulations. The analysis in the previous section on the relationships between macroscopic conductivities and microscopic flow

and thermal situations can still be applied here. Please note the change in interfacial conductance *C* does not affect the flow field and thus we can focus on the temperature response to different *C* values in Figure 7. As conductance *C* decreases, the solid domains become more insulated, resulting in a nearly constant temperature in each solid patch (i.e., the same colors and no isotherm contours in solid domains). With the solid part being insulated and thus its high-conductivity influence reduced, the fluid temperature around the solid domains exhibits a faster change along the thermal gradient direction. For example, the fluid temperatures near the solid surfaces at the narrowest gap location in Figure 7a are <sup>−</sup>0.0471 (left) vs. 0.0473 (right) for *<sup>C</sup>* <sup>=</sup> 10, <sup>−</sup>0.0364 vs. 0.0366 for *<sup>C</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> , and <sup>−</sup>0.0251 vs. 0.0254 for *<sup>C</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> . As discussed in the previous section, the large temperature difference at *C* = 10 introduces a larger tortuosity conductivity *Ktor*, and vise versa, as shown in Figure 8a. On the other hand, the dispersion conductivity *Kdis* increases as we reduce the interfacial conductance *C*, and the *C* influence becomes negligible for large aspect ratios. This change should be attributed to the temperature change in circulation areas; however, due to the complexity in flow structures and temperature changes with *C* in Figure 7, we are not able to provide a direct qualitative explanation here. Overall, the influence of the interface conductance *C* on the apparent conductivities is less dramatic in Figure 8 than that from the unit aspect ratio *β* in Figure 6, but it cannot be neglected for accurate analysis.

**Figure 7.** Simulation results of the flow streamlines (first column) and temperature distributions for different interfacial conductances *<sup>C</sup>* <sup>=</sup> 10 (second column), 5 <sup>×</sup> <sup>10</sup>−<sup>4</sup> (third column), and 5 <sup>×</sup> <sup>10</sup>−<sup>5</sup> (last column) with the aspect ratios *β* = 0.44 ((**a**) in the top row) and 1.44 ((**b**) in the bottom row) at *Re* = 50.

#### *4.3. Effect of the Porosity ǫ*

Unlike the effects of aspect ratio *β* and the interfacial conductance *C*, extensive studies have been conducted for the influence of the porosity *ǫ* on the macroscopic thermal behaviors of porous materials [21,61,62]. As mentioned in Section 1, previous pore-scale simulations usually set isothermal conditions on a pair of unit boundaries to create the thermal gradient and applied adiabatic conditions on other unit side surfaces [24,27]. Such artificial, unphysical treatments interfere with the natural heat transfer processes among microscopic periodic units, and the results from such boundary methods could be inaccurate or misleading [28]. In this section, we use the generalized periodic condition Equation (7) to examine the porosity effect on macroscopic conductivities *Ktor* and *Kdis*. Our next simulations cover a range of porosity of *ǫ* = 0.5∼0.99 for aspect ratios *β* = 0.69 and 1, and *ǫ* = 0.65∼0.99 for and *β* = 1.44. The Reynolds number is set at *Re* = 50 and the interfacial conductance is fixed at *C* = 10. Further reducing *ǫ* will increase the solid cylinder diameter *D* and thus make the gap between the solid surfaces very narrow (Figure 9). As discussed in Section 4.1, such a narrow gap may turn the flow unsteady for

the large aspect ratio *β* = 1.44, and also a finer spacial resolution is required to accurately capture the flow and thermal variations in the gap regions.

**Figure 8.** Plots of *Ktor* (**a**) and *Kdis* (**b**) changing with the aspect ratio *β* at three interfacial conductance values: *<sup>C</sup>* = 10 (black squares), 5 <sup>×</sup> <sup>10</sup>−<sup>4</sup> (blue circles), and 5 <sup>×</sup> <sup>10</sup>−<sup>5</sup> (red triangles). The Reynolds number is *Re* = 50.

**Figure 9.** Simulation results of the flow streamlines (first and third rows) and temperature distributions (second and last rows) for different porosity values (labeled on top of each column) with interfacial conductance *C* = 10 and Reynolds number *Re* = 50.

Results from this set of simulations are collected in Figure 9 for the flow and temperature distributions and Figure 10 for the calculated conductivities *Ktor* and *Kdis* at different porosity values. Please note that these figures are plotted in terms of the solid fraction 1 − *ǫ* as in Ref. [21]. With the solid diameter *D* increase (from left to right in Figure 9), a pair of circulation vortices are developed in the wake region above the solid cylinder, and the vortex size grows gradually till it completely fills the vertical space between two cylinders. Meanwhile, the original relatively organized temperature field is gradually distorted by the increasing size of the solid domain. Such changes in the flow and thermal patterns therefore affect the macroscopic thermal behaviors as characterized by the tortuosity (*Ktor*) and dispersion (*Kdis*) conductivities (Figure 10). The continuous increase in *Ktor* with the solid fraction 1 − *ǫ* is mainly because of the larger area (length in our 2D system) of the solid–fluid interface Γ (see Equation (14)). This trend is similar to that

observed in Ref. [21]; however, the *Ktor* magnitude is smaller here. This can be attributed to the different solid domain shapes: For the same porosity, compared to the circular shape in this work, the square solid domain shape in Ref. [21] yields a longer interface length, and more profoundly, half of the interface has the local normal direction aligned the *x* direction. All these are favorable for a larger *Ktor* according to its definition in Equation (14). As for the dispersion conductivity *Kdis* in Figure 10b, in general, *Kdis* increases with the solid fraction 1 − *ǫ*; however, unlike the monotonic growth in Ref. [21], local maximum and minimum states are observed here. Due to the complexity in flow and temperature fields as well as their involvement in the *Kdis* calculation, it is difficult to provide detailed insights and mechanism for the *Kdis* variations. Nevertheless, this could be an interesting topic to explore in the future. The general increasing trend can be qualitatively related to the increasing size of the circulation region, which is the major contributor to *Kdis* via the large *x* velocity component in this region. A similar analysis can be applied for the gentle variation and slow recovery of *Kdis* for 1 − *ǫ* = 0.2∼0.5 at *β* = 0.69 in Figure 10b, since the flow pattern remains almost unchanged in these systems (see the subplots for 1 − *ǫ* = 0.2 and 0.35 at *β* = 0.69 in Figure 9). Finally, for the three curves with different aspect ratios, we see *Ktor* is smaller but *Kdis* is larger for a higher aspect ratio *β*. This agrees well with our findings and discussions for the aspect ratio effect on conductivity in Section 4.1 (see Figure 6).

**Figure 10.** Plots of *Ktor* (**a**) and *Kdis* (**b**) changing with the solid fraction 1 − *ǫ* at three aspect ratio values: *β* = 0.69 (black squares), 1 (blue circles), and 1.44 (red triangles). The Reynolds number is *Re* = 50.

#### **5. Summary**

In this study, we have first justified and validated the generalized periodic condition [21,28] for pore-scale simulations of thermal flows in porous media. Extensive simulations have then been carried out using a 2D porous model to investigate the influence of several microscopic parameters on macroscopic thermal conductivity. Among them, the effects of the aspect ratio of the periodic unit and the thermal conductance at the pore surface have not been addressed adequately in previous studies. Our results show that these microscopic properties can dramatically change the flow and thermal fields in the microscopic porous structure, and affect the apparent thermal performances of the porous materials at the macroscopic level. Therefore, these microscopic factors need to be considered carefully for more accurate and reliable simulation results, which are crucial for both fundamental research and practical applications. In addition, thorough discussions are attempted to qualitatively explore the relationship between the apparent conductivity at the macroscopic level and the complex thermal flow situations in the microscopic porous

structure, and our analysis and findings could be helpful for a better understanding of the underlying thermal processes.

We are aware that several serious limitations exist for this study and one should not over interpret the results obtained from a specific system. The simple 2D porous model and the relatively low Reynolds numbers considered here may appear less realistic for practical systems; however, they are helpful for understanding the micro–macro relations and fundamental mechanisms involved in the complex thermal flow processes in porous materials. We have fixed the solid and fluid substance properties (conductivity, diffusivity, and the Prandtl number), which can certainly affect the apparent thermal behaviors as well. The 2D unit geometry adopted in this work is also symmetric along both the flow and thermal gradient directions; and the anisotropic effects could be an interesting topic for future research. The periodic conditions Equations (8) and (9) utilized in these calculations require the flow and thermal fields to be fully developed and thus they are not applicable to the entrance and development regions [63,64]. Moreover, in our simulations, we have assumed that the material properties are constant in the flow and heat transfer processes and thus steady or quasi-steady states can be established. In some situations, the microscopic porous structure may undergo a dynamic change due to swelling and erosion [65]; and caution must be taken for applying results from this study to such systems. Nevertheless, the boundary method, simulation results, and analysis discussions can be useful for the research and applications of porous thermal flows. The generalized periodic boundary condition, although presented in 2D, can be readily applied to three-dimensional pore-scale thermal flow simulations.

**Author Contributions:** Conceptualization, M.J. and J.Z.; methodology, M.J. and J.Z.; data curation, M.J.; formal analysis, M.J. and J.Z.; writing, M.J. and J.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Natural Science and Engineering Research Council of Canada (NSERC).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** M.J. acknowledges the financial support from the Queen Elizabeth II Graduate Scholarship in Science and Technology (QEII-GSST) at Laurentian University. The calculations have been enabled by the use of computing resources provided by Compute Canada (www.computecanada.ca).

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

#### **References**


#### *Article*
