Numerical Investigation on Two-Phase Flow Heat Transfer Performance and Instability with Discrete Heat Sources in Parallel Channels

With the rapid development of integrated circuit technology, the heat flux of electronic chips has been sharply improved. Therefore, heat dissipation becomes the key technology for the safety and reliability of the electronic equipment. In addition, the electronic chips are distributed discretely and used periodically in most applications. Based these problems, the characteristics of the heat transfer performance of flow boiling in parallel channels with discrete heat source distribution are investigated by a VOF model. Meanwhile, the two-phase flow instability in parallel channels with discrete heat source distribution is analyzed based on a one-dimensional homogeneous model. The results indicate that the two-phase flow pattern in discrete heat source distribution is more complicated than that in continuous heat source distribution. It is necessary to optimize the relative position of the discrete heat sources, which will affect the heat transfer performance. In addition, compared with the continuous heat source, the flow stability of discrete heat sources is better with higher and lower inlet subcooling. With a constant sum of heating power, the greater the heating power near the outlet, the better the flow stability.


Introduction
With the rapid development of integrated circuit technology, the heat flux of electronic chips has been sharply improved [1]. Therefore, heat dissipation has become the vital technology for electronic equipment and a more efficient radiator is required. Compared with single-phase flow, two-phase flow heat transfer capability is greater and the temperature of the heat transfer surface is more uniform, which makes flow boiling one of the most promising cooling technologies [2].
As a result, numerous investigations on heat transfer characteristics of flow boiling have been conducted in the past few decades. Yen et al. [3] adopted an experimental method to investigate the flow boiling phenomenon in microchannels with different cross sections. They indicated that the microchannel with square cross section has the best heat transfer capability because the square cross-section can provide a more effective nucleus of boiling. Deng et al. [4] compared the heat transfer performance of flow boiling between the microchannels with unique Ω-shaped cross sections and the conventional rectangular microchannels with the same hydraulic diameter by experiment. They found that the microchannels with unique Ω-shaped cross-sections can enhance the heat transfer because the distribution of liquid film is more uniform in annular flow and the narrow slot on the upside of the microchannels can suppress nucleated bubbly flow and early transform to the annular flow, which also contribute to the enhancement of two-phase heat transfer characteristics. Tran et al. [5] developed a formula to calculate the heat transfer coefficient based on a series of experimental data. The new correlation is applicable for smooth tubes with hydraulic diameters of 2.40-2.92 mm for the R113 (heat flux q = 8.8-90.8 kW/m 2 ; mass flux G = 50-400 kg/m 2 s), R12 (heat flux q = 7.5-129 kW/m 2 ; mass flux G = 44-832 kg/m 2 s), and R134a (heat flux q = 2.2-49.5 kW/m 2 ; mass flux G = 33-502 kg/m 2 s). Meng et al. [6] studied the flow boiling of R141b in vertical tubes and serpentine tubes through experiment and numerical simulation, and the results showed that the heat transfer can be enhanced by inclining the pipe. Huang et al. [7] studied numerically the flow pattern and heat transfer characteristics of flow boiling in the minichannels. They pointed out that the smaller the pipe diameter, the better the heat transfer. Catherine et al. [8] used a reduced-order model to simulate the flow boiling of R245fa and compared it experimentally. They focused on the effect of the time relaxation constant and r on flow and heat transfer, where they defined r as a function of the streamwise location. Daniel et al. [9] analyzed flow boiling in a microgap with circular pin fins using the Lee model and found that the evolution of two-phase flow is accelerated with the time relaxation constant increasing, which causes more uniform distribution of temperature. Similarly, Lee et al. [10] used the Lee model for analyzing thermal performance of the microcooler module in GaN HEMTs. The results showed that the heat transfer coefficient, temperature, and pressure drop depend on different mesh sizes and the time relaxation constant. Konstantinos et al. [11] adopted a numerical method to study the effect of surface wettability on flow boiling characteristics within microchannels. The results showed indeed that surface wettability plays a significant role on the HTC, and the hydrophilic and hydrophobic conditions are about 43.9% and 17.8% higher than the single-phase reference simulation, respectively. Li et al. [12] simulated the bubble departure characteristics and found that the surface tension is vital to simulate bubble departure diameter and active core density. Ali et al. [13] studied the subcooled flow boiling of different concentrations of alumina nanoparticles in a microchannel heat sink. Zhuan et al. [14] investigated the nucleate boiling of water in microchannels by using the VOF model. They displayed the whole process of bubbles changing in the microchannel, such as generating, growing, departing, and so on. Liu et al. [15] discussed the bubble growth and merger of refrigerant in a microchannel by using the coupled level set and volume of fluid method (CLSVOF). Yu et al. [16] investigated heat transfer and pressure drop of mixture refrigerant in a vertical rectangular minichannel, pointing out that as for heat transfer and pressure drop, the influence of heat flux was slight. Yang et al. [17] studied flow boiling of R141b in horizontal coiled tube, based on the Lee model, and the results showed that the temperature was significantly affected by the phase distribution. Using a numerical method, Steffen et al. [18] designed a channel in which vapor bubbles grew only along a predefined direction. Based numerical simulation results, Ma et al. [19] proposed a method to estimate the mean diameter of dispersed phase in saturated boiling.
However, in addition to the single heat source heating, there are also cases of multiple heat source heating used in the thermal engineering of equipment such as multichip electronics. Cho et al. [20,21] investigated the heat transfer effect of microchannel heat exchangers with variable cross-section when heated by nine discrete chips, and they found that the microchannel with a trapezoidal joint and progressively wider cross-section has the best heat transfer effect. Bhowmik et al. [22] used circulating cooling water for thermal management of four-chip electronic equipment and they emphasized that the relationship between convective heat transfer coefficient and Reynolds number is quite different from that of traditional single heat source heating. Kurose et al. and Islam et al. [23][24][25] experimentally studied the heat transfer characteristics and the flow distribution of refrigerant in nonuniform heating parallel channels. In order to improve the temperature uniformity of multiheating components, Ghaedamini et al. [26] proposed and designed a Y-shaped bifurcated microchannel network, and claimed that the smaller the Reynolds number, the more uniform the temperature distribution.
Although two-phase flow heat transfer has the advantages of strong heat transfer ability and uniformity, it has the trouble of flow instability. Two-phase flow instability will lead to flow oscillation, which will cause a series of damages such as mechanical vibration and heat transfer deterioration. Many experimental and theoretical studies on density wave oscillation (DWO), which is one of the most general instabilities, have been carried out. Colombo et al. [27] investigated the influence of bypass and channel inclination on DWO using RELAP5/MOD3.3. The results suggested that a sufficiently large bypass is equivalent to a constant pressure drop boundary, where the system stability is independent of bypass cross-sectional area. In addition, the flow stability is enhanced with increasing channel inclination. Based on the frequency domain method and two-phase model including void fraction, the numerical results of moving node scheme (MNS) and fixed node scheme (FNS) calculation of DWO were evaluated (Paruya et al. [28]). In their studies, MNS has higher computational efficiency and better convergence compared with FNS, and Froude number (Fr) causes different nonlinear bifurcations. Papini et al. [29] numerically analyzed DWO with the frequency domain method, and compared the calculation results of RELAP5 using a homogeneous flow model and COMSOL using the drift-flux model. It can be concluded that the frictional pressure drop of two-phase flow reduces the system stability, and appropriate bypass can effectively provide constant pressure drop. Lu et al. [30] researched DWO in parallel channels under axial nonuniform heating using the time domain method and homogeneous flow model, which showed that top-peaked heat flux can enhance the system stability in the whole region of inlet subcooling. Furthermore, the effects of system pressure and inlet resistance coefficient on flow stability under cosine heat flux were studied. A nodalized reduced order model (NROM) was developed to analyze DWO in forced circulation and natural circulation (Paul et al. [31]), which found that the stability boundaries of both circulations correlate with the number of nodes, and the Friedel model for friction factor calculation produces the most accurate result compared with the other four models. Paul et al. [32] studied the influence of system pressure on DWO in single channel with the frequency domain method and homogeneous flow model, and found that the generalized Hopf (GH) point, which divides subcritical and supercritical Hopf regions, disappears with system pressure, emphasizing that there is only the subcritical region in the stability boundary.
It can be seen from these studies that there are few existing studies on two-phase flow with discrete heat sources, although the heat transfer characteristics and flow instability of two-phase flow is widely studied. In this paper, in order to investigate the cooling of multielectronic chips, the heat transfer characteristics of flow boiling in parallel channels with discrete heat sources are studied numerically. Based on the one-dimensional homogeneous model, the two-phase flow instability with discrete heat sources heating is analyzed by programming. The results of this study are expected to provide some optimization methods for the application of two-phase flow cooling.

Physical Model
Based on literature research, it can be seen that the two-phase-flow heat transfer effect of the square channel is better, so this paper designed a square parallel channel to dissipate four discrete heat sources. As shown in Figure 1, the parallel channels formed in the aluminum plate are symmetrically distributed along the x, y, and z axes in space. The central axis of the branch channel is parallel to the central axis of the discrete heat source. In addition, the four discrete heat sources are symmetrically distributed along the x axis. The geometric parameters are shown in Table 1.

Parameter Value
Volume of aluminum plate (c × a × b)/cm 3 14 × 7 × 1 Volume of discrete heat sources (c2 × a2 × b2)/cm 3 2 × 2 × 0.3 Volume of single continuous heat source/cm 3 14 × 7 × 0.3 Cross section area of inlet (a1 × b1)/cm 2 1 × 0.5 Cross section area of outlet (a1 × b1)/cm 2 1 × 0.5 Cross section area of branch (a3 × a3)/cm 2 0.5 × 0.5 Axis distance of parallel channel (d1)/cm 4 The heating power of each discrete heat source is 60 W, and the solid walls are cooled by natural convection of air expected the solid walls of extended section, which are adiabatic. The ambient temperature is 15 • C. The inlet boundary is the mass flow boundary, which has a total mass flow of 0.015 kg/s, and the subcooling degree is 1.5 • C. The outlet boundary is the pressure boundary and operating pressure is 8 bar. Refrigerant R134a was used in the simulation; the properties of R134a are presented in Table 2. The heat source can be treated as heating copper.
where α, ρ, and → v are phase fraction, density, and velocity; S l and S v are the volumetric mass source terms; S ml and S mv are the mass transfer between the liquid phase and the vapor phase, which can be obtained in the phase-change model. Subscript "l" and "v" mean liquid phase and vapor phase.
Momentum conservation equation [6][7][8][9][10]13,17]: where P, µ, and t are pressure, viscosity, and time; → g and → F represent gravitational acceleration and interface-induced volume force. The surface tension of the interface is calculated by the continuum surface force model (CSF) model proposed by Brackbill et al. [33] According to the divergence theorem, the force at the surface can be expressed as a volume force used in the momentum conservation equation. Thus, → F is written as follows: where σ lv is surface tension coefficient, and κ v and κ l is surface curvature, which is defined as follows: Energy conservation equation: where E, k, and T are energy per kilogram, coefficient of thermal conductivity, and temperature, and S e is the volumetric energy source term.

Phase Change Model
The phase change model is proposed to provide the mass transfer S ml and S mv and the energy source term S e between the liquid phase and the vapor phase due to the boiling. In this paper, the evaporation-condensation phase change model proposed by Lee [34] was used.
According to the theory of gas dynamics, the mass flux through the phase interface during the phase transition can be expressed as follows: where β is the correction factor because the experimental theoretical evaporation is different. M is the molecular weight and R is the universal gas constant. T sat means saturation temperature. In order to maintain equilibrium between the liquid and gas phases, the chemical potential energy should be equal, so the equation between pressure and temperature can be obtained using the Clausius-Clapeyron equation: where h lv is latent heat and P sat means saturation pressure. Assuming that the bubbles have the same diameter d p , the area density of the phase interface is given as: Therefore, the mass transfer between the liquid phase and the vapor phase can be written as: To simplify the equations, the time relaxation factor parameters λ l and λ v are defined as: The mass transfer between the liquid phase and the vapor phase, S ml and S mv , can be written as follows: Thus, the energy source term can be calculated by the following equation: The time relaxation factors λ l and λ v should have appropriate values. If the value is too large, the calculation diverges easily and if the value is too small, the interface temperature deviates too much from the saturation temperature. Actually, the time relaxation factor is usually set as 0.1-100, referring to the studies by Lee [34] and Yang [17]. In this paper, λ l and λ v was set as 100 s −1 in order to maintain numerically the interface temperature within ±1 K compared with saturation temperature.

Solution Methods
Based on the pressure-based segregated solver with PISO scheme, the finite element method (FEM) was applied for solving the governing equations. The Green-Gauss nodebased method was chosen to solve the gradient algorithm PRESTO!, which was selected for pressure dispersion. In order to ensure accuracy, the momentum and energy equations are discretized by the second-order upwind scheme. When the residual of all equations is less than 10 −4 , the simulation is considered as convergence and the implicit body force is opened to ensure the convergence of numerical simulation.

Validation of the Simulated Method
In order to validate the accuracy of the simulation model, the same numerical model (VOF and evaporation-condensation phase change model) and same solution method (PISO) used in the present study were adopted to simulate the experimental study by Lin et al. [35]. Lin investigated two-phase heat transfer to R141b in the minichannel heated at the heat flux of 30,000 W/m 2 . According to the work of Lin, the boundary conditions of inlet and outlet were velocity boundary (0.42 m/s) and pressure boundary, and the number of cells is 160,979. Based on simulation and experimental data, the variation of heat transfer coefficient along the channel is shown in Figure 2, and the definition of the heat transfer coefficient is written as follows: where q, T wall , and T are the heat flux, average temperature of the wall, and the average temperature of fluid in the cross-section. From Figure 2, it can be seen that the deviation between experimental results and numerical value is less 15%, which shows that the simulated method is fully validating.
The grid sensitivity analysis on the result is indicated in Table 3. In order to ensure the solution accuracy and save on calculation time simultaneously, Mesh2 was selected.

Results and Discussion
In order to explore the influence of the relative position of discrete heat sources on the heat transfer characteristics of two-phase flow, the relative distances on the z-axis for two discrete heat sources in the same branch channel were changed to 3, 5, 7, 9 and 11 cm (d = 3/5/7/9/11 cm, respectively) and controlled other factors unchanged. Meanwhile, the group in which the single continuous heat source heated the aluminum plate with the same total heating power was treated as the control group. The average temperature of the single continuous heat source is 331.76 K. It can be seen from Figure 3 that the average temperature of four discrete heat sources is higher than that of the continuous heat source. The reason why the heat transfer effect of the cold plate is better when heated by a continuous heat source is because the heat exchange area of the continuous heat source is larger than that of the discrete heat source, and the refrigerant is heated more evenly. In addition, Figure 3 also shows that the heat transfer effect of two-phase flow is different when the relative position of the discrete heat source is different. In this paper, when the z-axis relative distance of two discrete heat sources in the same branch channel is 7-9 cm, the heat transfer effect of the two-phase flow is the best, and the heat transfer effect is the worst when the discrete heat source distribution is closest. When the z-axis relative distance of two discrete heat sources in the same branch channel is 7-9 cm, the two heat sources are close to the inlet and outlet of the channel, where the flow pulsation is more intense because the fluid direction changes and the fluid converges and disperses. Furthermore, it can be seen from Figure 4 that it is easier to nucleate and form bubbles at the corner of the flow channel. Therefore, the heat transfer capability of the inlet and outlet of the parallel channel is strong. In addition, the relative distance between discrete heat sources is large, so that the total heat will not be too concentrated, resulting in heat transfer deterioration. Thus, when the z-axis relative distance of two discrete heat sources is 7-9 cm, the temperature of the heat source is the lowest. On the contrary, when the distribution of discrete heat sources is closer, the more concentrated the total heat. If the heat is too concentrated, then it is easy to produce numerous bubbles and to form a gas film, which is not conducive to heat dissipation. When cooling the multiple chips, the temperature uniformity is also an important evaluation index. Uneven temperature distribution will generate thermal stress and affect the service life of the equipment. In Figure 5, "Heat source h1, h2, h3, h4" along the x-axis means four discrete heat sources; the location of the four discrete heat sources can be seen in Figure 1. The ordinate of Figure 5 represents the average temperature of heat sources. As shown in Figure 5, the two-phase flow heat transfer is helpful to improve the temperature nonuniformity and when the z-axis relative distance of two discrete heat sources in the same branch channel is 9 cm, the temperature distribution of the heat sources is most uniform.

Mathematical Model
Numerical methods for studying flow instability can be divided into two types: frequency domain method and time domain method. It should be noted that the frequency domain method can significantly shorten the time for numerical calculation through Laplace transform and linearization of the equation, but it cannot solve complex nonlinear problems. In this paper, the time domain method, referring to the investigations of Lu et al. [18] and Ma et al. [36], is adopted to analyze the flow instability.
Compared with heat transfer between fluid and channels, the investigation on flow instability pays greater attention to the axial changes of parameters instead of the one on cross-section. In order to effectively simplify the calculation, the following reasonable assumptions are made for the model:

1.
It is simplified to one-dimension, and only the variation of the parameters in the axial direction is considered; 2.
The entire tube is in thermal equilibrium; 3.
The fluid at the inlet of the pipeline is always in a state of supercooling; 4.
Subcooled boiling is not taken into account in the subcooling zone of the pipeline. The simplified physical model is shown in Figure 6.
Momentum conservation equation: where k is defined as: Energy conservation equation: When the fluid is in two-phase state: where ρ, u, h, p, x, and α are density, velocity, enthalpy, pressure, mass quality, and void fraction of fluid, respectively. The subscript g represents vapor, and the subscript l represents liquid. In addition, t means time and z refers to the axial coordinate, and q, g, λ, D e , k, and ∆p L are thermal power per unit volume, gravitational acceleration, Darcy friction coefficient, hydraulic diameter, the throttling coefficient at the inlet or outlet, and local resistance pressure drop, respectively. A staggered grid is used in this paper, in which p, ρ, and h are stored at the center of the control volume, while u is stored at the boundary of the control volume, as shown in Figure 7. The first order upwind difference and semi-implicit scheme are employed to discretize the equation, and the following equations can be obtained by the chain rule of partial derivatives.
Mass conservation equation: Momentum conservation equation: Energy conservation equation: where the superscript n and the subscript i represent the next time level and the control volume sequence number, respectively, while ∆t and ∆z represent the time step and the space step.
Assuming that all parameters of the current time level are known, combining and simplifying Equations (24)- (27) gives the following form: where δp = p n − p is the pressure difference between the next and the current time level. In addition, A j (j = 1,2,3) and B are constants. In order to eliminate the step of flow distribution and improve the calculation speed in the program, the upper header and the subheader are divided into two control volumes, and then the discretization is carried out in the same way as noted above. The matrix equation of the whole system can be obtained by integrating Equation (33) for every control volume: δp in Equation (32) can be calculated by the Gaussian elimination method, and then the remaining unknowns, including u n , h n , ρ n , can be determined by substituting δp back into Equations (24) and (28)-(30), respectively. The numerical computations noted are implemented by FORTRAN codes.

Solution Method and Model Validation
The transient calculation is performed while an initial field of system is given. After running for a sufficiently long time, the system reaches steady state when the flow rate of each tube does not change over time and the thermal disturbance is now applied, that is, the thermal power of a tube increases by 0.1%. After a short time, e.g., 0.01 s, the thermal power is restored to its original state, which is equal to a pulse consisting of two step functions, and then the changes of flow rate in the two tubes are observed.
If divergent oscillation occurs, that is, the amplitude of the oscillation of two tubes becomes increasingly larger, that means the thermal power is too high. Otherwise, it is a damped oscillation. Until the oscillation becomes a constant amplitude, the thermal power at this time is the critical thermal power. As shown in Figure 8, when the degree of inlet subcooling is 30 K and the total thermal power is 4730.6 W, the parallel channel oscillates in constant amplitude as a result of the thermal disturbance noted above according to the parameters set in Table 4.  By comparing the set model with experimental data [37], it can be found that when the pressure of the system is 140 bar and the inlet undercooling degree is 120 K, different total critical thermal power per unit area of system can be obtained for different total mass flow, as shown in Figure 9. The errors between the calculated values and the experimental values are within ±15%, which indicates that the two are in excellent agreement, and the correctness of the model is verified.

Result and Discussion
In the study of flow instability, the two-dimensional plane composed of dimensionless numbers, N pch and N sub , is generally used to characterize the stability boundary. N pch means the phase change number and N sub represents the subcooling number. Their definitions follow: where Q, W, and h fg refer to heating power, mass flow rate, and latent heat of vaporization. Additionally, ∆h in = h l − h in is undercooling enthalpy of the inlet. The physical properties and geometric parameters as shown in Table 3 are used for calculation, and the flowing medium is water. First, the influence of discrete and continuous heat source distribution on flow instability under the same total heating power is analyzed. The continuous distribution of the heat source is shown in Figure 6 while the setting of the discrete heat source for comparison is shown in Figure 10. In order to be representative, the discrete heat sources are set to be uniformly distributed along the pipe segment, and each pipe has two discrete heat sources. Calculated according to the parameters in Table 3, the stability boundary finally obtained is shown in Figure 11. When the inlet undercooling degree is relatively low and relatively high, the stability boundary of the discrete heat source is located on the right of the continuous heat source boundary, which indicates that the flow of discrete heating is more stable than that of continuous heating. However, when the inlet undercooling is moderate, the stability of the continuous heat source is better than that of the discrete heat source. The analysis shows that when the inlet undercooling degree is low, the fluid is heated by the discrete heat source with high heat flux and directly turns into single-phase steam. When the degree of undercooling at the inlet is high, the enthalpy of the fluid is still at a low value even after passing through the discrete heat source, and the nonheated section of the channel is basically in the undercooled region. Regardless of whether the fluid in most channels is in gas or liquid phase, its density does not change dramatically, thus inhibiting the occurrence of flow instability. However, when the degree of undercooling at the inlet is within a specific range, resulting in the fluid heated by the discrete heat source just being in the two-phase region, the nonheated region is equivalent to extending the two-phase region, that is, the region with periodic changes of density is prolonged, causing the deterioration of the flow stability.
The numerical distribution of discrete heat sources also affects the flow instability. For the discrete heat source shown in Figure 10, keeping the position of the heat source and the total heating power of the single tube unchanged, but changing the respective proportion of total heating. power The working conditions are set as shown in Table 5. Heat source 1 and heat source 2 refer to the discrete heat source near the inlet and outlet, respectively. A stability boundary diagram of three working conditions is shown in Figure  12. The results show that the numerical distribution of the discrete heat source, which can affect the flow stability, is not only related to the uniformity but also relevant to the distribution sequence along the flow path. In the stability boundary diagram, the order from left to right is working conditions 2, 1, and 3, which is also the order in which the stability gradually improves. It can be considered that the greater the power of the discrete heat source close to the inlet and the smaller the power close to the outlet, the worse the flow stability will be, and vice versa. This can be attributed to the fact that the high-power discrete heat source close to the inlet rapidly increases the enthalpy of the fluid, resulting in drastic changes in the density of the fluid in most of the subsequent flow channels, thus promoting the occurrence of flow instability. However, when the high-power heat source is near the outlet, the enthalpy value of most fluids in the flow passage is low and the density changes little. Even if the density changes dramatically after heating by the heat source, the remaining flow passage is not enough to fully show the flow instability. Thus, flow stability is improved when the high-power discrete heat source is closer to the flow outlet or the low-power discrete heat source is closer to the flow inlet.

Conclusions
In this paper, the two-phase heat transfer characteristics and instability in parallel channels with discrete heat source arrangement were studied by numerical simulation. The main conclusions are as follows: (1) The relative positions of discrete heat sources affect the heat transfer effect of twophase flow, and the heat transfer performance can be improved by reasonably designing the relative position of heat sources. (2) The closer the distribution of discrete heat sources, the worse the heat transfer effect; for the working condition designed in this paper, the heat transfer effect is best when the distance between the centers of the two discrete heat sources on the same branch channel are within a range of 7-9 cm. (3) Compared with the continuous heat source, the discrete heat source with uniform distribution can enhance the flow stability under low and high inlet undercooling. (4) When the high-power discrete heat source is closer to the flow outlet and the lowpower discrete heat source is closer to the flow inlet, the flow stability is improved. Funding: This research received no external funding.

Institutional Review Board Statement:
The study did not involve humans or animals.

Informed Consent Statement:
The study did not involve humans.

Data Availability Statement:
The study did not report any data.

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