Simulation Study on the Non-Uniform Characteristics of Boiling Flow and Heat Transfer in Parallel Small Channels

1  Introduction

Currently, thermal management has evolved from a supportive measure ensuring normal equipment operation into a core bottleneck restricting further improvements in system performance. The average heat flux of certain advanced chips has reached the kilowatt-per-square-centimeter level, with even higher values in local hot spots. Meanwhile, global data center electricity consumption has exceeded 1% of total power consumption, with approximately 40% dedicated to cooling systems [1]. Severe heating issues lead to degraded device performance, deteriorated reliability, and shortened lifespan. They have also become a thermal barrier constraining the development of strategic industries, including high-performance computing, artificial intelligence, power electronics, 5G/6G communications, and radar. Traditional single-phase cooling technologies, limited by their heat transfer capacity, struggle to meet dissipation demands. In contrast, flow boiling heat transfer, due to its high efficiency in utilizing latent heat, is regarded as a highly promising next-generation thermal management solution. Therefore, there is an urgent need to conduct research on flow boiling within parallel small channels.

Compared to the simple model of a single small channel, parallel channels offer higher heat dissipation flux. By increasing the number of channels, the heat dissipation area can be expanded to meet the cooling requirements of high-power devices. However, their flow stability is highly sensitive and susceptible to inlet-outlet pressure coupling, leading to system-level flow and pressure oscillations [2]. Furthermore, their heat transfer characteristics are dominated by flow distribution. The uneven distribution results in significant variations in channel quality, flow patterns, and Critical Heat Flux (CHF) among channels, causing local overheating. Lee and Mudawar [3] comprehensively summarized non-uniformity in parallel small channels boiling, CHF enhancement mechanisms, and types of flow instability. They explicitly pointed out that the multi-channel coupling effect was the core feature distinguishing them from a single channel. Deng et al. [4] confirmed through visualization that in multi-channels, bubble coalescence paths are interfered with by vapor from adjacent channels. This interference leads to the formation of vapor bridges. In contrast, bubble motion in a single channel remains independent and free from such inter-channel interference. This was the microscopic origin of heat transfer non-uniformity. Xiong et al. [5] emphasized that parallel flow passages were widely used in heat exchangers, evaporators, and boiling or condensation systems. The capacity of the system throughput can be increased, and the heat and mass transfer rates can be enhanced. Vontas et al. [6], through numerical studies, compared a single rectangular microchannel with a hydraulic diameter of 200 μm and a parallel microchannel heat sink with four channels. In the parallel configuration, each channel has a hydraulic diameter of 50 μm. They found that the parallel configuration reduced average thermal resistance by 23.3%, whereas the single channel achieved only a 5.4% reduction. This quantitative comparison fully reveals the significant thermal performance advantages of multi-channel systems. Karayiannis and Mahmoud [7] discussed in detail experimental studies on flow boiling in single tubes and rectangular multi-channels. They explicitly stated that flow instability, backflow, and their impact on heat transfer rates are fundamental problems unique to multi-channel systems. Zhang et al. [8] clearly pointed out that uneven two-phase flow distribution in parallel channels was typically associated with uneven flow distribution within the inlet manifold and instabilities in uniform distribution. Aka and Narayan [9] were the first to apply entropy balance analysis to predict two-phase flow distribution in parallel channels. They discovered that as the flow rate decreased, uniform distribution became severely non-uniform, accompanied by a significant change in entropy generation rate exceeding 30%. While entropy generation analysis for single channels considers only their own irreversibility, the entropy balance for multi-channels reveals a different behavior. The system tends to select flow distribution patterns that minimize entropy generation. Additionally, Yang et al. [10] and Tuo and Hrnjak [11] argued that flow reversal of the working fluid within parallel microchannels caused pressure drop oscillations, resulting in non-uniform flow distribution.

In recent years, numerous reports have been published regarding the non-uniform distribution of flow and heat transfer in flow boiling within parallel multi-channels. Liu et al. [12] experimentally investigated gas-liquid two-phase slug flow distribution in horizontal parallel microchannels. The results revealed the intrinsic link between slug flow characteristics (such as slug length and velocity) and uneven flow distribution. Kurose et al. [13] conducted experiments on two parallel small channels with non-uniform heating. They found that when the heat flux difference between channels was large, significant flow maldistribution occurred. The channel with higher heat flux received less than the average flow rate, and its heat transfer coefficient (HTC) near the outlet decreased. Miglani et al. [14] experimentally studied the impact of non-uniform heating on flow boiling distribution in parallel microchannels. The results indicated that heat flux differences between channels significantly exacerbated flow maldistribution. Yang et al. [15] experimentally investigated the flow boiling heat transfer characteristics of deionized water in rectangular microchannels with large aspect ratios. They found that the HTC gradually decreased as channel size, heat flux, saturation pressure, and mass flow rate decreased. Moreover, the enhancing effect of increased mass flow rate on HTC diminished as channel size decreased. Al-Zaidi et al. [16,17] experimentally studied the flow boiling heat transfer and two-phase flow characteristics of HFE-7100 in parallel microchannels. They discovered that HTC peaks when vapor quality approached zero and subsequently decreased with increasing vapor quality. They also noted that HTC was barely affected by mass flow rate but was significantly influenced by heat flux. Furthermore, through Fast Fourier Transform (FFT) analysis, they concluded that the frequency and amplitude of pressure drop fluctuations were positively correlated with heat flux. Theoretical modeling of complex thermo-fluid phenomena provides essential insights for experimental and numerical studies. Abdelsalam et al. [18] developed a semi-analytical solution for Magnetohydrodynamics and heat transfer of a couple-stress fluid subject to thermal radiation in a porous medium, using the homotopy perturbation method. Their findings contribute to understanding how multiple physical fields, including magnetic, thermal, and porous effects, interact in fluid systems.

Ledinegg instability [19] is the most classic theory explaining the generation of non-uniform flow distribution in flow boiling within parallel multi-channels. When the characteristic curve between pressure drop and flow rate of a channel exhibits a negative slope, the uniform flow state becomes unstable. Consequently, the system spontaneously transitions to a stable state with uneven flow. At the inlet manifold of parallel channels, gas-liquid two-phase flow is prone to separation due to differences in density and inertia, leading to vastly different void fractions at the inlets of individual channels. This constitutes the origin of non-uniformity. Subsequently, the thermal-fluid-solid coupling effects [20] exacerbate this non-uniformity. The thermal coupling between channels (via conduction through a shared substrate) significantly influences the disparity. Meanwhile, minor differences in heating power can be drastically amplified by flow instability. In systems with multiple parallel channels, Ledinegg instability may lead to non-uniform flow distribution. As a result, some channels receive less flow compared to uniform distribution conditions, experience severe dryout, trigger premature CHF, and compromise the reliability of the heat sink. Several measures have been proposed to mitigate two-phase flow maldistribution, including inlet restrictors and throttling valves [21], increasing system pressure [22], and active flow control mechanisms at channel inlets [23]. Barzegar Gerdroodbary et al. [24] numerically investigated the convective heat transfer of ferrofluid in a serpentine tube. They found that magnetic fields can induce secondary vortices, which improved the fluid mixing.

Currently, research on flow boiling and heat transfer non-uniformity in parallel small channels still faces several key limitations. These limitations severely restrict the engineering application and theoretical development of this technology. Existing studies mostly rely on experimental methods to analyze non-uniform flow and heat transfer characteristics within small channel systems. However, due to limitations in measurement techniques, obtaining comprehensive information on flow and thermal fields remains difficult. Furthermore, current descriptions of non-uniform phenomena in flow and heat transfer distribution across channels are relatively isolated. The degree of both non-uniformities has not been reasonably quantified. There is a lack of clear understanding regarding the coupling relationship between them. Therefore, this paper employs numerical simulation methods to investigate the flow boiling and heat transfer characteristics of R134a. The focus is specifically on the non-uniform behavior of these processes within parallel, rectangular small channels. The differential characteristics of flow and heat transfer parameters among channels under varying mass flux and heat flux conditions were systematically investigated. The quantitative characterization methods for the degree of flow and heat transfer non-uniformity were proposed. Ultimately, the predictive correlations for the non-uniformity degrees of flow and heat transfer were constructed. This study can provide theoretical tools for the quantitative assessment of flow and heat transfer non-uniformity in multi-channel systems.

2  Methodology

2.1 Grid Geometry and Boundary Conditions

This study employs the Finite Volume Method (FVM) to investigate the non-uniform characteristics of flow boiling and heat transfer in parallel small channels. Fig. 1 illustrates the schematic of the computational domain and boundary conditions. The model dimensions of the test section are 20 mm (length) × 10 mm (width) × 0.5 mm (height), comprising flow channels and fins, with a total of 15 channels. The light yellow regions represent the flow channels, each measuring 20 mm × 0.5 mm × 0.5 mm, while the gray regions denote the fins with a width of 0.2 mm. The bottom surface, highlighted in red, serves as the heating wall. Its geometric parameters align with the test model described in experimental literature [25]. To eliminate inlet effects and ensure fully developed flow upon entering the test section, an entrance section of 10 mm in length is added at the inlet. Similarly, an exit section of 10 mm is incorporated at the outlet to mitigate backflow effects, minimizing their impact on the flow field within the test section. A solid boundary is extended in the negative z-direction from the leftmost channel to simulate heat exchange between the actual channel and the solid side section in the simulation. Since the simulation results show that the averaged temperature of the heating wall remains almost unchanged when the width of the side section is greater than 5 mm, its width is determined as 5 mm. In the simulation setup, the entrance section, exit section and flow channels are defined as fluid domains, whereas the fins and side section are designated as solid domains. Given the presence of 15 channels and the symmetric nature of the entire configuration, a symmetry plane (40 mm × 0.5 mm) is established at the mid-width. This approach reduces the grid number and enhances computational efficiency without compromising simulation accuracy. Consequently, this study analyzes only 7 channels, numbered 1 through 7 along the positive z-axis. The boundary conditions applied in this study are illustrated in Fig. 1. A velocity inlet boundary condition is specified at the inlet of the entrance section. The inlet velocity is calculated based on the mass flux at three levels of 400, 600, and 800 kg/(m2·s). The inlet temperature is fixed at 20°C, which corresponds to a subcooled liquid state at the operating pressure of 770 kPa. The bottom surface of the flow channels, highlighted in red in Fig. 1. It is subjected to a constant heat flux boundary condition with four heat flux levels ranging from 200 to 800 kW/m2 applied uniformly over the entire heated surface. The side solid section and the top surfaces of the channels and fins are treated as adiabatic walls since heat transfer from these surfaces is negligible compared to the bottom heating wall. A symmetry boundary condition is applied at the mid-width plane of the domain. This reduces the computational cost by half, as the geometry and flow are symmetric about this plane. Finally, all solid-fluid interfaces are set as no-slip walls with coupled thermal boundary conditions, ensuring continuity of heat flux and temperature across the interface. Detailed numerical values are provided in Tables 1 and 2.

Figure 1: Schematic diagram of calculation domain and boundary conditions.

Fig. 2 illustrates the grid conditions of the calculation domain. A structured grid is employed for the spatial discretization of the domain. Since the entrance and exit sections involve simple flow and the side section only involves heat conduction, a coarser grid resolution is sufficient for accurate simulation of the test section. In contrast, the flow channels and fins, which constitute the heated test section, are discretized with a finer grid to enhance the calculation accuracy of the gas-liquid two-phase flow. Furthermore, as shown in the magnified view, the fluid region adjacent to the heating wall undergoes additional local refinement. This strategy allows for a more precise resolution of critical phenomena such as the ONB, bubble nucleation, and boundary layer heat transfer. Following this grid refinement, the dimensionless wall distance (y+) near the heated surface is maintained below 0.02.

Figure 2: Grid conditions of the calculation domain.

2.2 Fluid and Solid Properties

In this study, R134a is selected as the working fluid, and aluminum is chosen as the solid wall material. The thermophysical properties of both the liquid and vapor phases of R134a, as well as those of aluminum, are assumed to be constant. The primary thermophysical parameters for the gas, liquid, and solid phases within the computational domain, obtained from the REFPROP database, are listed in Table 3.

2.3 Numerical Model

This section describes the numerical models used in the simulation. The simulation process consists of four main steps. First, the VOF model tracks the gas-liquid interface by solving the continuity equations for each phase. Second, the CSF method models surface tension at the interface. Third, the energy equation is solved to obtain the temperature field. Fourth, the Lee model simulates the phase change process when the local temperature reaches the saturation temperature. The SST k-ω turbulence model is also employed to capture turbulence effects in the boiling flow.

2.3.1 Governing Equations

In this study, the Volume of Fluid (VOF) model is employed to track the free interface between different phases in the boiling flow. The spatiotemporal evolution of the volume fraction for each phase is described by solving the continuity equations for the two-phase volume fractions. Specifically, the liquid and vapor phases are designated as the primary and secondary phases, respectively. The corresponding continuity equations are presented as Eqs. (1) and (2).

∂(αlρl)∂t+∇⋅(αlρlu→)=Sl(1)

∂(αvρv)∂t+∇⋅(αvρvu→)=Sv(2)

The subscripts l and v denote parameters for the liquid and vapor phases, respectively. Sl is the mass transfer rate from liquid to vapor due to phase change, Sv is the mass transfer rate from liquid to vapor due to phase change. The sum of αl and αv within each control volume is always equal to unity.

The Continuous Surface Force (CSF) method is employed to model the surface tension effects at the curved gas-liquid interface during the flow boiling process. The surface tension is converted into a volumetric force and incorporated as a source term into the momentum equation, as shown in Eq. (3).

∂(ρu→)∂t+∇⋅(ρu→u→)=−∇p+∇⋅μ[∇u→+(∇u→)T]+ρg→+F→vol(3)

where the F→vol is the volumetric force due to surface tension, computed according to Eq. (4).

F→vol=σαlρlκl∇αl+αvρvκv∇αv0.5(ρl+ρv)(4)

where σ is the surface tension coefficient. κl and κv denote the liquid and vapor curvatures. Their definitions are given by Eq. (5).

κl=−κv=Δαl|∇αl|(5)

The ∇αl is calculated using the Least-Squares Curvature Reconstruction method. This approach enhances numerical robustness for geometries with high aspect ratios. Furthermore, an interface compression scheme is applied, coupled with grid refinement in near-wall regions and around channel corners, to mitigate curvature estimation errors induced by steep gradients.

The continuous liquid and vapor phases share a common energy equation, as presented in Eq. (6).

∂(ρE)∂t+∇⋅[u→(ρE+p)]=∇⋅(λ∇T)+SE(6)

where E represents the mass-averaged internal energy. It is defined in Eq. (7). λ is the thermal conductivity. SE is the energy source term due to phase change.

E=αlρlEl+αvρvEvαlρl+αvρv(7)

Under conditions with significant temperature differences between the two phases, the calculation accuracy of the interfacial temperature field is severely compromised. Large temperature gradients introduce anisotropic coefficients into the governing equations, leading to substantial numerical errors. Therefore, a volume-fraction-weighted interpolation method is employed to determine the thermophysical properties ϕ of the fluid, including thermal conductivity λ, dynamic viscosity μ, and density ρ, as expressed in Eq. (8).

ϕ=αlϕl+αvϕv(8)

The thermal behavior of the solid domain is governed by the transient heat conduction equation, as shown in Eq. (9).

ρcp∂T∂t=∇⋅(λ∇T)(9)

where cp is the specific heat capacity of the solid.

2.3.2 Turbulence Equation

The SST k-ω turbulence model is employed to simulate boiling turbulence within parallel small channels. This model accounts for low-Reynolds-number effects and combines the advantages of both the k-ω and k-ε models, demonstrating unique effectiveness in capturing turbulence characteristics and vortex dissipation phenomena during phase change processes [26].

2.3.3 Phase-Change Model

The Lee model is employed to simulate the phase change process between the liquid and vapor phases. This model utilizes the temperature difference between each computational cell and the saturation temperature as the primary driving force for phase change. Specifically, liquid evaporation occurs when the local liquid temperature exceeds the saturation temperature (T ≥ Tsat), whereas vapor condensation occurs otherwise. The mass transfer rate between phases is governed by Eq. (10) and is incorporated as a mass source term into the continuity equations given in Eqs. (1) and (2).

Sv=−Sl={rlαlρlT−TsatTsat,T≥TsatrvαvρvT−TsatTsat,T<Tsat(10)

where rl and rv represent the evaporation and condensation coefficients, respectively. They serve as critical tuning parameters to regulate model accuracy and significantly influence the simulation results of the phase change process [27]. Increasing the r value may lead to numerical divergence, whereas decreasing it can result in an excessive temperature difference between the phase interface and the saturation temperature. Regarding the selection of r, it is generally based on empirical experience from previous studies and specific case requirements. It is widely accepted that the default r value is appropriate when the temperature difference between the phase interface and the saturation temperature does not exceed 1 K [28]. In this study, by comparing the numerical simulation results of boiling flow with experimental data, the values of rl and rv were ultimately set to 50 and 2000, respectively.

The phase change process involves heat absorption and release. The corresponding energy source term, SE, is determined by Eq. (11).

SE=Slhfg=−Svhfg(11)

where hfg is the latent heat of vaporization.

With the above models, the governing equations for mass, momentum, energy, and turbulence are solved iteratively at each time step. The specific solver settings are detailed in the following section.

2.4 Solution Methods

To simplify the computation and conserve computational resources, this study first performs a steady-state calculation for the flow within the channels. Once all flow field data stabilize, the simulation switches to a transient calculation. The specific solver settings are detailed in Table 4. In the transient solution, the time term is discretized using the First-order implicit scheme, which offers unconditional stability and is well-suited for the intense phase change processes characteristic of boiling flow. Gradient reconstruction employs the Least Squares Cell Based method, enabling accurate capture of local gradient variations at the gas-liquid interface. Pressure-velocity coupling is achieved via the Coupled algorithm. This fully implicit approach effectively avoids the pressure-velocity decoupling issues commonly encountered in transient two-phase flow simulations. The pressure term is discretized using the Pressure Staggering Option (PRESTO!) scheme to accurately calculate pressure gradients near the phase interface. The momentum and energy equations utilize the Bounded Central Differencing scheme, which maintains second-order accuracy while employing limiters to prevent non-physical oscillations. Turbulence parameters (turbulent kinetic energy and specific dissipation rate) are also discretized using the Least Squares Cell Based method. To ensure the accuracy of phase change rate calculations, the convergence criterion for residual is set to 10−6. The simulation proceeds to the next time step when the residual meets this criterion or when the maximum number of iterations (20 steps) is reached.

2.5 Validation of Grid and Time Step Independence

Numerical computation time and accuracy are directly correlated with the grid number and the time step. While increasing the grid number and decreasing the time step generally improve calculation accuracy, they also significantly increase computational cost. Fig. 3 illustrates the impact of grid number and time step on the average temperature Tw of the heating wall under operating conditions of a mass flux G = 400 kg/(m2·s) and a heat flux q = 400 kW/m2. The black and red solid lines represent the results of the grid and time step independence tests, respectively. The grid independence study was performed using eight grid levels ranging from 0.4 to 7.2 million cells. The average temperature of the heating wall was chosen as the monitoring parameter because it is sensitive to near-wall resolution. As shown in Fig. 3, the temperature variation between successive grid refinements decreased progressively. The relative change between 4.0 and 4.98 million cells was 0.24%, and between 4.98 and 7.8 million cells, it was less than 0.1%. Based on this analysis, the grid with 4.98 million cells was selected for all subsequent simulations.

Figure 3: Average temperature of the heating wall calculated with different grid numbers and time steps.

Regarding the time step independence test, simulations were performed with the 4.983-million-cell across four time steps. As shown in Fig. 3, as the time step decreases exponentially, the variation in the average temperature of the heating wall diminishes. When the time step is reduced below 10−4 s, the average temperature remains virtually unchanged. Therefore, a time step of 10−4 s was selected. The chosen grid number and time step are highlighted by dashed boxes in the figure.

2.6 Validation of Numerical Model

In order to validate the reliability of the numerical model, both qualitative and quantitative comparisons were conducted between the simulation results and the experimental data from Ref. [25]. Fig. 4 compares the numerical simulation results of gas-liquid two-phase flow patterns in parallel small channels with visualization images under a mass flux condition of G = 400 kg/(m2·s). The simulation perspective in Fig. 4a is aligned with the field of view captured by the high-speed camera in Fig. 4b to ensure consistency in the comparison. Under different heat fluxes, bubbly flow, slug flow, and churn flow were observed in the parallel multi-channels. In the bubbly flow regime, bubbles are distributed discretely, with continuous coalescence and attachment of small bubbles near the wall. In the slug flow regime, bubbles within the channels grow larger. A liquid slug that contains numerous smaller bubbles appears, connecting large bubbles constrained by the channel geometry. As the length of the liquid bridges continuously shortens and the vapor slugs elongate, significant bubble breakup occurs, leading to the churn flow regime characterized by fragmented bubbles. Given the high degree of consistency between the simulated flow patterns and the experimental visualization results, the feasibility and accuracy of the proposed model are qualitatively validated.

Figure 4: Comparison between simulated gas-liquid flow patterns and corresponding visualization images under G = 400 kg/(m2·s).

Fig. 5 compares the average HTC of the heating wall obtained from numerical simulations and experiments under conditions of heat fluxes q = 200–800 kW/m2 and mass flux G = 400 kg/(m2·s). The error band of ±30% is indicated by red dashed lines. The simulated average HTC of the heating wall was calculated using Eq. (12). The experimental data are reproduced from Ref. [25] for validation purposes. As shown in the figure, all simulated HTC values fall within ±30% across the entire heat flux range. This error level is considered acceptable for predictive studies of multiphase flow, as similar or even higher error ranges have been reported in the literature [29,30]. The main sources of error in this study include the assumption of constant thermophysical properties and the neglect of heat conduction in the solid wall, which can be addressed in future work by incorporating temperature-dependent properties and conjugate heat transfer modeling. At q = 200–300 kW/m2, the simulation slightly overpredicts HTC. The slight overprediction at low heat fluxes may be due to the Lee model’s sensitivity to subcooled boiling initiation. At q = 600–800 kW/m2, a slight underprediction is observed. The underprediction at high heat fluxes may be attributed to the simplified bubble coalescence and breakup models that do not fully capture the complex churn flow dynamics. Despite these minor discrepancies, the overall agreement is satisfactory. This quantitative comparison confirms that the numerical model is capable of predicting flow boiling heat transfer in parallel small channels with reasonable accuracy.

hsim=qTw−Tl(12)

where Tl represents the fluid temperature at the first grid layer adjacent to the heating wall.

Figure 5: The comparison between the simulated and experimentally measured HTCs of heating wall under mass flux of G = 400 kg/(m2·s).

Several limitations of this study should be acknowledged. The numerical model is based on the VOF-Lee approach. This model relies on empirical coefficients. These coefficients were calibrated under specific conditions. They may not be universally applicable. The current investigation is limited to a specific channel geometry. The channel dimensions are fixed at 0.5 mm by 0.5 mm with a rectangular shape. The operating parameters are also restricted. The mass flux ranges from 400 to 800 kg/(m2·s). The heat flux ranges from 200 to 800 kW/m2.

3  Results and Discussion

3.1 Non-Uniform Flow Characteristics

Fig. 6 shows the influence of heat flux on the gas-liquid two-phase flow pattern and wall temperature distribution in parallel small channels under the condition of mass flux G = 400, 600 and 800 kg/(m2·s). The simulation results with low mass flux G = 400 kg/(m2·s) are taken as representatives to analyze the bubble flow state and evolution characteristics. As shown in Fig. 6a, a large number of bubbles appear in each channel at all heat fluxes, and the flow pattern in each channel remains basically the same. When q = 200 kW/m2, the bubbles present discrete spherical and ellipsoidal shapes and move downstream close to the heating wall. Along the flow direction, the bubble size gradually increases, but there is no obvious coalescence, which conforms to the slug flow morphology. At the same time, the temperature distribution on the heating wall is uniform, and the channel demonstrates superior overall thermal performance. When the heat flux increases to q = 400 kW/m2, discrete bubbles grow rapidly in the upstream of the channel, presenting a slug flow pattern. As the fluid flows downstream and continues to be heated, multiple large bubbles approach, collide and merge with each other. These merged bubbles continue to move and grow in the channel until their size is limited by the channel height and width, and then they begin to be axially elongated in the flow channel. At this time, the elongated bubbles continue to extend longitudinally along the channel under the drive of buoyancy and flow, which conform to the typical Taylor bubble structure and present a slug flow pattern. At the outlet section of the channel, the slug bubbles break, and a local high-temperature area appears on the bottom heating wall. When the heat flux increases to q = 600 kW/m2, the slug flow region in the channel is further shortened, and the breaking position of the slug bubbles moves upstream. As a result, the area of the local high-temperature area expands, and the difference among each channel also increases. At the highest heat flux q = 800 kW/m2, the gas-liquid interface completely loses stability, and there are no clearly distinguishable slug bubbles. The gas-liquid two-phase presents the boiling and churning flow, conforming to the churn flow pattern. The local high-temperature area covered by the churn flow bubbles at the outlet section of the channel further expands.

Figure 6: Liquid-gas flow patterns and wall temperature distributions in multi-channels under different heat fluxes and mass fluxes.

Fig. 6b,c shows the flow and heat transfer in multi-channels under the higher mass fluxes. As the mass flux increases, the bubble nucleation position gradually moves downstream. There are fewer obvious discrete bubbles in the channel because the discrete bubbles coalesce on the heating wall and present a thin strip-shaped bubble morphology in the flow direction under the action of the mainstream inertial force. When the thin strip-shaped bubbles gradually become thicker in the flow direction and approach the channel size, slug bubbles appear in the channel and form a slug flow pattern.

Due to differences in the flow pattern of each channel, it is necessary to analyze the flow rate differences of each channel. Fig. 7a–c shows the flow rate distribution of different channels under the condition of mass flux G = 400, 600, 800 kg/(m2·s). The dashed lines in the figures are the ideal average mass flow rates of each channel under the corresponding mass fluxes. The error bars are the ranges of fluctuations in the actual average mass flow rates. It can be seen from the figures that there are fluctuations and uneven distributions of the flow rate among the channels. Under the same mass flux, as the heat flux increases, the fluctuations and uneven distributions of the flow rate also increase. This occurs because higher heat flux intensifies bubble nucleation and growth. It leads to a larger gas void fraction in channels with slightly higher local temperatures. The presence of vapor locally increases the two-phase flow resistance in those channels. This positive feedback mechanism amplifies the initial small differences among channels, resulting in more severe flow maldistribution. As the mass flux increases, the flow rate distribution among the channels gradually becomes more uniform, and the influence of the heat flux on the flow rate distribution also weakens. This is because higher mass flux enhances the inertial force of the incoming liquid. The increased inertia suppresses the relative impact of localized vapor-induced resistance variations. Consequently, even when bubble generation is intense, the stronger inertia helps maintain a more balanced flow distribution across channels, making the system more resilient to heat flux disturbances. In Fig. 7a, the ideal average mass flow rate of each channel is about 0.133 g/s, and the fluctuation range of most channels is within 0.12–0.14 g/s. However, the difference between Channel 1 and other channels is relatively large. At q = 400 and 800 kW/m2, the mass flow rate of the channel deviates significantly from the average value. As the heat flux increases, the degree of deviation of each channel also increases. In Fig. 7b, the ideal average mass flow rate of each channel is 0.2 g/s, and the fluctuation range of most channels is within 0.19–0.21 g/s. Similarly, the difference between Channel 1 and other channels is relatively large, but the flow rate distribution of each channel is more concentrated than that in Fig. 7a. The fluctuation range is also significantly reduced, indicating that as the mass flux increases, the distribution tends to be more balanced, the stability of the system is improved, and the disturbing effect of the heat flux is suppressed. In Fig. 7c, the ideal average mass flow rate of each channel is 0.267 g/s. The fluctuation range of all channels is within the range of 0.26–0.28 g/s. The data under the four heat fluxes basically coincide, indicating that the influence of the heat flux on the flow rate distribution is very small for the largest mass flux condition. Even under high heat fluxes, the data does not show obvious deviation.

Figure 7: Mass flow fluctuation in each channel, increase with heat flux, decrease with mass flux.

Fig. 8a–d shows the variation of the cross-sectional gas void fraction of each channel with the axial position under different heat fluxes at the mass flux G = 400 kg/(m2·s). It is difficult to determine the critical position of the flow pattern transition through the flow pattern in Fig. 6. Therefore, it is necessary to quantitatively distinguish the dominant regions of each flow pattern. The previous research [25] showed that for the R134a working fluid, when the gas void fraction was 0 < αv < 0.35, the flow pattern belonged to the bubbly flow. When it was 0.35 < αv < 0.65, the flow pattern belonged to the slug flow. When it was 0.65 < αv < 0.9, the flow pattern belonged to the churn flow. For reference, the demarcation dotted lines of αv = 0.35 and 0.65 are drawn to better find the critical position of the flow pattern transition. Overall, as the heat flux increases, the rising rate of the gas void fraction of each channel accelerates. Meanwhile, the fluctuations of the curves also increase significantly. This behavior is caused by the intermittent nature of two-phase flow at higher heat fluxes. As bubbles grow and coalesce into larger vapor slugs, the instantaneous gas void fraction oscillates between liquid-rich and vapor-rich states. The amplitude of these oscillations increases as the flow regime transitions. The flow evolves from bubbly flow with stable and low fluctuations to slug flow with periodic and moderate fluctuations, and finally to churn flow with chaotic and high fluctuations. The increasing fluctuation magnitude reflects the growing instability of the gas-liquid interface. In Fig. 8a, the overall gas void fraction is within 0.35, and the flow pattern is the bubbly flow, which is consistent with the flow pattern map. In Fig. 8b, the gas void fraction of most channels breaks through 0.35 near x/L = 0.6 and enters the slug flow. The gas void fraction of Channel 1 is always lower than that of the other channels. This is because the mass flow rate of Channel 1 is the largest relative to the other channels under this condition, as can be seen from Fig. 7a. In Fig. 8c, the gas void fraction of multiple channels exceeds 0.35 near x/L = 0.4 and quickly enters the slug flow. The gas void fraction of Channel 1 is significantly higher than that of the other channels and enters the churn flow near x/L = 0.5. However, the gas void fraction of the other channels only slightly exceeds 0.65 at the end of the channel and enters the churn flow. The mass flow rate of Channel 1 is the smallest relative to the other channels in Fig. 7a, so its gas void fraction is relatively high. In Fig. 8d, the gas void fraction of almost all channels has exceeded 0.35 near x/L = 0.3 and exceeds 0.65 near x/L = 0.6. Although the fluctuations of the curves increase significantly, the gas void fraction in each channel at the outlet is basically around 0.9.

Figure 8: Increase regularity of cross-sectional gas void fraction αv in axial direction for each channel under different heat fluxes and G = 400 kg/(m2·s).

Fig. 9a–d shows the variation of the cross-sectional gas void fraction for each channel with the axial position under different heat fluxes at a mass flux of G = 600 kg/(m2·s). Generally, as the heat flux increases, the rising rate of the gas void fraction in each channel speeds up. However, when the mass flux is increased, at the same heat flux, the gas void fraction in each channel decreases significantly at the same position, and the difference in the gas void fraction among channels as well as the fluctuations of each curve decrease. In Fig. 9a, the overall growth of the gas void fraction is slow. The gas void fraction remains around 0.2 at the outlet and the flow regime is bubble flow. In Fig. 9b, the gas void fraction in Channel 1 is always lower than that in other channels. The flow regime remains bubble flow at the outlet, finally. The gas void fraction in other channels reaches 0.35 around x/L = 0.85, and the flow regime transforms into slug flow. In Fig. 9c, the gas void fraction in Channel 1 is always lower than that in the other channels. This is because the mass flow rate in Channel 1 is the largest among the other channels as shown in Fig. 7b. The gas void fraction in the other channels reaches 0.35 around x/L = 0.65 and enters slug flow, but they do not reach 0.65 finally. In Fig. 9d, the gas void fraction in Channel 1 is always higher than that in the other channels, and only Channel 1 shows the churn flow regime. Because the mass flow rate in Channel 1 is the smallest among the other channels in Fig. 7c, the gas void fractions in the other channels are basically the same, reach 0.35 around x/L = 0.6. Once the slug flow occurs, the gas void fractions increase slowly and reach about 0.55 at the outlet.

Figure 9: Increase regularity of cross-sectional gas void fraction αv in axial direction for each channel under different heat fluxes and G = 600 kg/(m2·s).

Fig. 10a–d shows the variation of the gas void fraction in each channel with the axial position under different heat fluxes at a mass flux of G = 800 kg/(m2·s). Generally, with the increase in the heat flux, the rising rate of the gas void fraction in each channel accelerates, but their values are basically the same at the same position. At a high mass flux, compared with Figs. 8 and 9, the increase in the gas void fraction in each channel is slow, and there is no channel showing churn flow, which is consistent with the gas-liquid flow shown in Fig. 6c. In Fig. 10a,b, the gas void fraction is always lower than 0.35 and the flow pattern remains bubbly flow. The fluctuations become larger for the higher heat flux condition of q = 400 kW/m2. In Fig. 10c, the gas void fractions in each channel are basically the same, and slug flow only appears at the tail of the channel. In Fig. 10d, the gas void fraction in each channel remains basically unchanged in the x/L = 0–0.35 section, and then fluctuates. The gas void fraction reaches 0.35 at x/L = 0.65 and slug flow appears. It reaches about 0.5 at the outlet. Therefore, increasing the mass flux can effectively inhibit excessive local vaporization, delay the deterioration of the flow pattern, and enhance the ability of the system to resist high heat flux disturbances.

Figure 10: Increase regularity of cross-sectional gas void fraction αv in axial direction for each channel under different heat fluxes and G = 800 kg/(m2·s).

According to the axial distribution of gas void fraction under different working conditions, the axial regions dominated by each flow regime are plotted in Fig. 11. The blue, orange and red parts represent the bubbly flow, slug flow and churn flow. Overall, as the mass flux decreases and the heat flux increases, the position of the flow regime transition advances, and more flow regime types appear in the channel. Under the mass flux of G = 600 and 800 kg/(m2·s), almost all channels are in bubbly flow or initially transition to slug flow along the entire length, and there is no churn flow. The uniformity of flow regime distribution is higher than the mass flux of G = 400 kg/(m2·s). In Fig. 11a, when the heat flux is q = 200 kW/m2, the whole channel is occupied with bubbly flow. When the heat flux is increased to q = 400 kW/m2, slug flow appears near x/L = 0.65 in most channels, but the dominant flow regime in the channel is still bubbly flow. When the heat flux is q = 600 kW/m2, slug flow appears near x/L = 0.4 in most channels, and churn flow appears near x/L = 0.75. The dominant flow regime in the channel has changed to slug flow. Most of the channel is occupied by slug bubbles in Fig. 6a, presenting a typical slug flow regime. When the heat flux is q = 800 kW/m2, slug flow appears near x/L = 0.3 in the channel, and churn flow appears near x/L = 0.6. The dominant flow regime in the channel has changed to churn flow. The gas-liquid interface completely loses stability, and there are no clearly distinguishable slug bubbles, presenting a churn flow regime in Fig. 6a. In Fig. 11b, except for the Channel 1 under the highest heat flux condition of q = 800 kW/m2, where the dominant flow regime is slug flow, the churn flow appears near x/L = 0.85, the dominant flow regime in all channels and conditions is bubbly flow. In Fig. 11c, with the further increase of mass flux, only when the heat flux is q = 600 and 800 kW/m2, the slug flow appears at the tail of the channel, and the dominant flow regime in the channel is bubbly flow.

Figure 11: Axial distribution of gas-liquid flow patterns for each channel under different heat fluxes and mass fluxes.

The gas-liquid phase change can cause changes in the flow resistance of each channel, which in turn affects the pressure drop of each channel. Since the churn flow no longer appears in the channels and the dominant flow pattern in all channels is bubbly flow at the highest mass flux of G = 800 kg/(m2·s), the pressure drop characteristics of the channel for the mass flux of G = 400 and 600 kg/(m2·s) under different heat fluxes are compared in Fig. 12. The left and right figures show the time-domain and frequency-domain distributions of the pressure drop, respectively. The frequency-domain data are obtained by performing the FFT method on the time-domain data. The pressure drop data within the time period of 0.1–0.3 s are selected for analysis, which contains a complete flow development cycle.

Figure 12: Time-domain (left) and frequency-domain (right) diagrams of pressure drop fluctuations in each channel under different heat fluxes and mass fluxes.

In terms of the time-domain response, as the heat flux increases, the pressure drop of the channel increases. The fluctuation characteristics across all channels remain largely consistent. Notably, Channel 1 exhibits a marginally lower pressure drop at a heat flux of q = 400 kW/m2, which is attributed to its relatively low gas void fraction, as illustrated in the Figs. 8 and 9. Besides, the periodic pressure pulsations occur in each channel, with the amplitude gradually increasing as the heat flux increases. Bubbles nucleate, grow, and detach periodically on the heating wall, causing periodic changes in local resistance. As the mass flux increases, the pressure drop also increases, but the amplitude of the fluctuation decreases, and the regularity of the waveform is enhanced. In Fig. 12a, when q = 200 kW/m2, the pressure drop of the channel is almost a smooth straight line with a small amplitude of fluctuation. This indicates that the flow is in an approximately single-phase or incipient boiling state at this time, with few and evenly distributed bubbles generated, and the system is highly stable as a whole. As the heat flux increases to q = 400 kW/m2, small-scale, quasi-periodic pressure pulsations emerge in each channel. These pulsations exhibit a slight increase in amplitude, indicating that nucleate boiling has become fully developed. Further increasing the heat flux to q = 600 kW/m2, the pressure drop fluctuations intensify significantly, with not only a marked increase in amplitude but also a distortion in the waveform shape. When q = 800 kW/m2, the pressure drop signal shows highly unsteady characteristics, with a chaotic waveform, frequent large jumps, steep drops, and long-time troughs. These features are usually related to the accumulation of vapor slugs, local dryout-rewetting cycles, or strong interactions between gas-liquid two-phase flows. Such behavior indicates that the system is approaching or has entered the flow instability region, with the risk of heat transfer deterioration and local burnout. In Fig. 12b, when q = 200 and 400 kW/m2, the characteristics of the pressure drop are basically the same as those at low mass fluxes. When q = 600 kW/m2, periodic sawtooth-like fluctuations can be observed. However, compared with low mass fluxes, the amplitude of vibration is generally reduced, and the waveform is smoother and more repetitive. When the heat flux is further increased to q = 800 kW/m2, the fluctuation also increases at high mass fluxes. However, a certain periodicity is still retained as a whole, and the main oscillation mode is clearly distinguishable. These characteristics indicate that the system is still in a controllable strong boiling region and has not entered the dangerous heat transfer deterioration stage.

In terms of the frequency-domain response, as the heat flux increases, the low-frequency oscillation becomes stronger. As the mass flux increases, there is basically no obvious change at low heat fluxes. At high heat fluxes, it shows lower frequencies and higher-energy peaks. When G = 400 kg/(m2·s), q = 200, 400 kW/m2 and G = 600 kg/(m2·s), q = 200, 400, 600 kW/m2, the dominant flow pattern is bubbly flow. The spectral curve is basically flat, and the energy is concentrated at extremely low frequency, with no obvious peak. When G = 400 kg/(m2·s), q = 600 kW/m2 and G = 600 kg/(m2·s), q = 800 kW/m2, the dominant flow pattern is slug flow. A recognizable main frequency peak begins to appear in the frequency spectrum within the range of 5–10 Hz, which is consistent with the typical bubble detachment frequency range. Secondary harmonic peaks begin to appear in the interval of 15–25 Hz, which may be due to bubble coalescence, slug movement, or the reflection and superposition of local pressure waves in the channel. When G = 400 kg/(m2·s), q = 800 kW/m2, the dominant flow pattern is churn flow. The main frequency is still concentrated in 5–10 Hz. Energy is further amplified, with richer higher-order harmonic components. The frequency spectrum broadens to above 30 Hz, showing typical broadband random oscillation characteristics.

3.2 Non-Uniform Heat Transfer Characteristics

The non-uniformity in heat transfer among parallel channels arises from the coupling between flow distribution and local boiling dynamics. Uneven flow distribution reduces liquid supply in some channels. This increases the void fraction and allows vapor to cover the heating wall, forming a thermal barrier that lowers the heat transfer coefficient. Boiling dynamics, including bubble nucleation, growth, and flow regime transition, alter the gas-liquid interface distribution and wall wetting conditions. These changes further amplify or mitigate heat transfer differences among channels. Since the non-uniform flow may lead to differences in the HTC of each channel, the non-uniform of heat transfer characteristics was analyzed. Fig. 13 presents scatter plots of the HTC for each channel under varying heat fluxes at mass fluxes of G = 400, 600 kg/(m2·s). The plots reveal the coupled effects of heat flux and mass flux on thermal performance. There are differences in the HTC of each channel, with no obvious monotonic characteristics from 1 to 7. The HTC initially maintains a high level before decreasing significantly as heat flux increases. This is because the gas void fraction in the channel increases as the heat flux increases, separating the liquid from the heating wall. The thermal conductivity of steam is much lower than that of liquid; the gas layer becomes the main thermal resistance. With the increase of mass flux, the HTC of each channel decreases overall under a specific heat flux, and the maximum difference between channels decreases from 20 to 10 kW/(m2·°C). For Fig. 13a, the HTC of channel 1 is higher than that of the other channels. This difference remains almost unchanged as the heat flux increases. For Fig. 13b, under the higher mass flux, the HTC in the channels is more uniform. Under the condition of q = 200 kW/m2, it stays around 50 kW/(m2·°C). When q = 800 kW/m2, the HTC of channel 1 is lower than that of the other channels. This is because only channel 1 exhibits churn flow, as shown in Fig. 11.

Figure 13: Average HTC on heating wall for each channel under different heat fluxes and mass fluxes.

Fig. 14 shows the distribution of the local HTC at different heat fluxes under the mass fluxes G = 400 and 600 kg/(m2·s) along different axial positions in the channel. With the increase in heat flux, the overall HTC shows a downward trend. At the same heat flux, the HTC at a high mass flux is relatively higher than that at a low mass flux. This is caused by the increase in mass flux, which enhances the convective heat transfer intensity. When G = 400 and 600 kg/(m2·s) and q = 200 kW/m2, the HTC of all channels gradually increases from 50–60 kW/(m2·°C) at the inlet to 125–150 kW/(m2·°C) near x/L = 0.5, and maintains this value until the outlet. This is because only bubbly flow appears under these conditions as shown in Fig. 11. For G = 600 kg/(m2·s), q = 400 kW/m2, the HTC increases monotonically along the axis, which is caused by slug flow in the channel. When G = 400 kg/(m2·s), q = 400 kW/m2 and G = 600 kg/(m2·s), q = 600 and 800 kW/m2, the HTC first increases and then decreases along the axial direction. This is due to the expansion of the slug flow range in the channels. Under the lower mass flux condition of q = 400 kW/m2, the peak positions in each channel differ. At the higher mass flux, the HTC distribution in each channel tends to become uniform. As the heat flux increases, the HTC peak decreases from 145 to 100 kW/(m2·°C), and its position moves upstream from x/L = 0.9 to 0.65. When G = 400 kg/(m2·s), q = 600 and 800 kW/m2, the HTC remains at a relatively low level in the axial direction. This is caused by the increased range of churn flow control within the channels. The HTC differences among the channels are small, and heat transfer deterioration occurs in all of them. Based on the analysis of Figs. 11 and 14, it can be observed that the HTC differences among channels are small under bubbly flow and churn flow conditions, but large under slug flow conditions. This is because under bubbly flow, the bubbles are small and uniformly distributed, leading to similar flow conditions across channels. Under churn flow, intense gas-liquid mixing averages out the differences among channels. Under slug flow, gas slugs randomly and intermittently block the channels, causing significant variations in flow resistance and wall cooling effects across different channels.

Figure 14: Axial distribution of HTC on the heating wall for each channel under different heat fluxes and mass fluxes.

Since three flow patterns coexist under the mass flux G = 400 kg/(m2·s) as can be seen from Fig. 11, the relationship between the gas void fraction and the HTC in the channel under this mass flux and different heat fluxes is plotted in Fig. 15. The dashed line represents the changing trend of the HTC. With the increase of the gas void fraction, the HTCs under four heat fluxes all show a trend of first increasing and then decreasing. With the increase in the heat flux, the maximum gas void fraction in the channel increases from 0.3 to 0.89, but the maximum HTC decreases from 170 to 50 kW/(m2·°C). The gas void fraction corresponding to each HTC peak increases from 0.18 to 0.65. For regions dominated by bubbly flow (αv < 0.35).

Figure 15: Relationship between cross-sectional gas void fraction αv and HTC on heating wall under different heat fluxes and G = 400 kg/(m2·s).

3.3 Relationship between Flow and Heat Transfer Non-Uniformity

From the above analysis, it can be seen that there are differences in the flow and heat transfer characteristics of each channel. Moreover, this non-uniformity changes with the working conditions. In this paper, the coefficients of variation from statistics are used to quantitatively characterize the dimensionless numbers β1 and β2, which describe the non-uniformity of flow rate and HTC, as shown in Eqs. (13) and (14).

β1=1n−1∑i=1n(mi−mn¯)21n∑i=1nmi(13)

β2=1n−1∑i=1n(hi−h¯)21n∑i=1nhi(14)

Fig. 16 shows the non-uniformity of flow rate and HTC under different working conditions. In Fig. 16a, the β1 increases with the increase of heat flux. At the lowest mass flux G = 400 kg/(m2·s), with the increase of heat flux, the increase amount of β1 keeps basically the same. For a specified heat flux, the β1 shows a non-linear relationship of first rapidly decreasing and then slowly decreasing with the increase of mass flux. This is because, as the mass flux increases, the flow inertia in each channel is enhanced, the impact of flow resistance caused by gas-liquid phase change is weakened, and uniformity is improved. As for the non-uniformity of heat transfer shown in Fig. 16b, the β2 decreases almost linearly with the increase of mass flux. Furthermore, it also decreases with the increase in heat flux. However, when the heat flux increases from 400 to 600 kg/(m2·s), the decreasing rate of β2 increases. This operating condition change corresponds to the transition from slug flow to slug flow. As the heat flux further increases, the β2 continues to decrease. This is because the HTC of slug flow and churn flow is at a lower level, resulting in reduced heat transfer non-uniformity between channels.

Figure 16: Non-uniformity of mass flow rate and HTC for multi-channels under different heat fluxes and mass fluxes.

A good quantitative relationship can be observed between flow regime transition and β1. This is because flow regime transition tends to occur more readily under conditions of increasing heat flux and decreasing mass flux. Similarly, β1 also increases with increasing heat flux and decreases with decreasing mass flux, showing a consistent trend with flow regime transition. Specifically, when only bubbly flow exists in the channels without other flow regimes, β1 is less than 0.035. When the flow regime transitions from bubbly flow to plug flow, β1 ranges between 0.035 and 0.14. When the flow regime further transitions from plug flow to churn flow, β1 exceeds 0.14.

Fig. 17 shows the relationship between the β1 and β2 under different working conditions. The dotted line represents the fitted trend under the mass fluxes G = 400, 600 and 800 kg/(m2·s). It can be seen that the β2 decreases as the β1 increases under a specified mass flux. With the increase of heat flux, the rising rate of β1 is much smaller than the falling rate of β2. In contrast, the β2 increases as the β1 increases under a specified heat flux. Under a given mass flux, β1 only varies within the range of 0.05, and the heat flux has little effect on the degree of flow non-uniformity. However, it has a significant impact on the degree of heat transfer non-uniformity, with the variation range reaching 0.14 at each mass flux. As the mass flux increases, the span of β1 variation under different heat fluxes increases from 0.02–0.07 to 0.12–0.16, indicating that the mass flux has a significant effect on the degree of flow non-uniformity. However, the span of β2 variation under different heat fluxes only increases from 0.01–0.14 to 0.04–0.16, indicating that the mass flux has little effect on the degree of HTC non-uniformity.

Figure 17: Relationship between non-uniformities of mass flow rate and HTC on heating wall for multi-channels under different heat fluxes and mass fluxes.

3.4 Predictive Correlation for Flow and Heat Transfer Non-Uniformity

Previous text found that non-uniformity is significantly related to G and q, so this study hypothesizes that the degree of flow and heat transfer non-uniformity is mainly affected by Bo, Re, and We. The boiling number Bo (=q/(G·hfg)) characterizes the ratio of heat flux to the kinetic energy of flow, which reflects the intensity of phase change and captures the effect of thermal input on non-uniformity. The Reynolds number Re (=u·L/υ) represents the ratio of inertial force to viscous force. It governs the flow regime and determines the flow distribution characteristics among parallel channels. The Weber number We (=ρ·u2·L/σ) accounts for the shear force acting on bubbles. It captures the role of surface tension in bubble deformation and gas-liquid interface stability. Together, these three dimensionless numbers collectively describe how heat input, flow inertia, and surface tension influence the non-uniformity of flow boiling in parallel small channels. We construct an explicit predictive correlation for the degree of flow and heat transfer non-uniformity between parallel multiple channels based on Bo, Re, and We. As shown in Eqs. (15) and (16).

β1=f(Bo,Re,We)(15)

β2=f(Bo,Re,We)(16)

This paper uses a power-law function form to fit the correlation. In this function, β is taken as the dependent variable. The expression composed of the remaining dimensionless numbers characterized by ψ serves as the independent variable. The specific formula is as shown in Eqs. (17) and (18). Where k is an empirical constant, and n, a, b, c are unknown exponents. Using the numerical values of β1 and β2 under different working conditions and the least squares method, we can obtain the unknown parameters in the equations.

β=k⋅ψn(17)

ψ=Boa⋅Reb⋅Wec(18)

After fitting the simulation data, it is obtained that, for β1, k = 0.95, n = 2, for β2, k = 1.04 × 10−4, n = 1.17. These can help us draw curves to detect the margin of error in simulated data. As shown in Fig. 18, the data in this paper all fall within the 20% error band, indicating the reliability of the Eqs. (19) and (20). It can effectively predict the non-uniformity of boiling flow and heat transfer under the conditions where the mass fluxes ranged from 400 to 800 kg/(m2·s) and the heat fluxes ranged from 200 to 800 kW/m2.

β1=0.95⋅Bo0.5⋅Re−0.3⋅We0.2(19)

β2=1.04×10−4⋅Bo−0.85⋅Re−0.2⋅We0.15(20)

Figure 18: Comparison of proposed non-dimensional correlation, based Bo, Re, We and β1, β2, with simulation data.

The dataset covers mass fluxes of 400, 600, and 800 kg/(m2·s), heat fluxes of 200, 400, 600, and 800 kW/m2, resulting in 12 distinct operating condition combinations. To avoid overfitting during the fitting process, a leave-one-out cross-validation method was employed. With a total of 12 data points, 11 points were used for fitting and the remaining 1 point for validation. This procedure was repeated 12 times, with each data point serving as the validation set once. The cross-validation results confirmed that the prediction errors remained within the acceptable ±20% range for all validation sets, indicating that overfitting did not occur and the correlations have good generalizability.

The exponents of Re in both β1 and β2 are negative, so increasing Re can reduce both non-uniformities. This can be achieved most directly by increasing the mass flux. The exponents of We are positive, so decreasing We can reduce both non-uniformities, which can be achieved by reducing the channel hydraulic diameter. The exponents of Bo have opposite signs: positive in β1 and negative in β2. Therefore, a lower Bo improves flow uniformity, while a higher Bo improves heat transfer uniformity. The choice should be made based on the design priority.

4  Conclusion

In this paper, a comprehensive numerical simulation of the boiling flow and heat transfer characteristics in parallel small channels under different operating conditions is carried out. The non-uniform characteristics of the flow and heat transfer are studied. The flow and heat transfer characteristics and distribution differences of each channel are analyzed. The dimensionless numbers β1 and β2, representing the non-uniformity degrees of flow and heat transfer, are proposed. Their correlations with known operating conditions are fitted. The main conclusions are summarized as follows:

(1)   The VOF coupled with the Lee model can well simulate the boiling flow in parallel small channels. The error of the averaged HTC on the heating wall obtained by simulation is controlled within ±30% compared with the experimental results.

(2)   The non-uniformity degree of flow rate between parallel channels weakens with the increase of the mass flux and intensifies with the increase of the heat flux. At low mass fluxes, the dominant flow pattern in the channel changes from bubbly to slug and churn flow with the increase of the heat flux. However, the dominant flow pattern is basically bubbly flow for higher mass flux conditions.

(3)   The non-uniformity degree of HTC between parallel channels weakens with the increase of the mass flux and heat flux. Under a fixed mass flux, HTC first increases then decreases with gas void fraction, and its peak value declines as heat flux rises due to flow regime transitions.

(4)   As the non-uniformity degree of flow rate β1 increases, the non-uniformity degree of HTC β2 decreases under a specified mass flux, while it increases under a specified heat flux. The predicted errors of the fitted correlations of β1 and β2 are both controlled within ±20%.

The proposed predictive correlations for the non-uniformity degrees are empirically derived based solely on the present dataset. Their extrapolation to other channel sizes or working fluids requires further validation. Future work should address these limitations by extending the parameter range and incorporating more sophisticated modeling strategies.

Acknowledgement: Not applicable.

Funding Statement: This research was funded by the National Natural Science Foundation of China [No. 52576046], the Open Foundation of State Key Laboratory of High-End Compressor and System Technology [No. SKL-YSJ202406], the National Natural Science Foundation of China [No. 52376037], and the Zhejiang Provincial Key Research and Development Project [No. 2024C01117].

Author Contributions: Conceptualization, Chi Zhong; methodology, Chi Zhong; validation, Chi Zhong and Xiao Wang; formal analysis, Chi Zhong; investigation, Chi Zhong; resources, Yang Liu and Linmin Li; data curation, Xiao Wang and Bo Ye; writing—original draft preparation, Chi Zhong; writing—review and editing, Yang Liu; supervision, Yang Liu and Linmin Li; project administration, Yang Liu; funding acquisition, Yang Liu. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data will be made available on request.

Ethics Approval: Not applicable. It does not involve human or animal subjects.

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

Nomenclature

Roman symbols

Greek symbols

Subscripts

References

1. Wu Z, Xiao W, He H, Wang W, Song B. Jet-enhanced manifold microchannels for cooling electronics up to a heat flux of 3000 W cm−2. Nat Electron. 2025;8(9):810–7. doi:10.1038/s41928-025-01449-4. [Google Scholar] [CrossRef]

2. Ghiaasiaan SM. Two-phase flow, boiling and condensation in conventional and miniature systems. Cambridge, UK: Cambridge University Press; 2008. [Google Scholar]

3. Lee J, Mudawar I. Two-phase flow in high-heat-flux micro-channel heat sink for refrigeration cooling applications: part I—pressure drop characteristics. Int J Heat Mass Transf. 2005;48(5):928–40. doi:10.1016/j.ijheatmasstransfer.2004.09.018. [Google Scholar] [CrossRef]

4. Deng C, Luo XP, Feng ZF, Zhang RD. Research on boiling heat transfer characteristics and visualization of refrigerant in rectangular microchannels. J Refrig. 2015;36(6):1–5. (In Chinese). [Google Scholar]

5. Xiong T, Liu G, Huang S, Yan G, Yu J. Two-phase flow distribution in parallel flow mini/micro-channel heat exchangers for refrigeration and heat pump systems: a comprehensive review. Appl Therm Eng. 2022;201(7830):117820. doi:10.1016/j.applthermaleng.2021.117820. [Google Scholar] [CrossRef]

6. Vontas K, Andredaki M, Georgoulas A, Miché N, Marengo M. The effect of hydraulic diameter on flow boiling within single rectangular microchannels and comparison of heat sink configuration of a single and multiple microchannels. Energies. 2021;14(20):6641. doi:10.3390/en14206641. [Google Scholar] [CrossRef]

7. Karayiannis TG, Mahmoud MM. Flow boiling in microchannels: fundamentals and applications. Appl Therm Eng. 2017;115(3):1372–97. doi:10.1016/j.applthermaleng.2016.08.063. [Google Scholar] [CrossRef]

8. Zhang L, Hu G, Bi XT. Two-phase flow in parallel channels: mal-distribution, hysteresis and mitigation strategies. Chem Eng Sci. 2022;247(1):117044. doi:10.1016/j.ces.2021.117044. [Google Scholar] [CrossRef]

9. Aka T, Narayan S. Using entropy balance to determine multiphase flow distribution in thermally decoupled parallel channels with shared inlet and outlet headers. Phys Fluids. 2024;36(5):054118. doi:10.1063/5.0207373. [Google Scholar] [CrossRef]

10. Yang KS, Jeng YR, Huang CM, Wang CC. Heat transfer and flow pattern characteristics for HFE-7100 within microchannel heat sinks. Heat Transf Eng. 2011;32(7–8):697–704. doi:10.1080/01457632.2010.509774. [Google Scholar] [CrossRef]

11. Tuo H, Hrnjak P. Periodical reverse flow and boiling fluctuations in a microchannel evaporator of an air-conditioning system. Int J Refrig. 2013;36(4):1263–75. doi:10.1016/j.ijrefrig.2013.02.001. [Google Scholar] [CrossRef]

12. Liu Y, Sun W, Wang S. Experimental investigation of two-phase slug flow distribution in horizontal multi-parallel micro-channels. Chem Eng Sci. 2017;158:267–76. doi:10.1016/j.ces.2016.10.021. [Google Scholar] [CrossRef]

13. Kurose K, Miyata K, Hamamoto Y, Mori H. Flow boiling heat transfer and flow distribution in non-uniformly heated parallel mini-channels. In: Proceedings of the 3rd Thermal and Fluids Engineering Conference (TFEC); 2018 Mar 4–7; Fort Lauderdale, FL, USA. p. 1499–502. doi:10.1615/tfec2018.che.021663. [Google Scholar] [CrossRef]

14. Miglani A, Soto A, Weibel JA, Garimella SV. The effect of uneven heating on the flow distribution between parallel microchannels undergoing boiling. J Electron Packag. 2021;143(4):041107. doi:10.1115/1.4052532. [Google Scholar] [CrossRef]

15. Yang XQ, Qiu SZ, Jia XH, Yin HF, Jia DN, Lu DH. Experimental study on the heat transfer of water boiling flow in narrow horizontal rectangular channels. Nucl Power Eng. 2007;28(3):38–42. (In Chinese). doi:10.1115/icone10-22212. [Google Scholar] [CrossRef]

16. Al-Zaidi AH, Mahmoud MM, Karayiannis TG. Flow boiling of HFE-7100 in microchannels: experimental study and comparison with correlations. Int J Heat Mass Transf. 2019;140(17–18):100–28. doi:10.1016/j.ijheatmasstransfer.2019.05.095. [Google Scholar] [CrossRef]

17. Mahmoud MM, Karayiannis TG. Pool boiling review: part II—heat transfer enhancement. Therm Sci Eng Prog. 2021;25(1):101023. doi:10.1016/j.tsep.2021.101023. [Google Scholar] [CrossRef]

18. Abdelsalam SI, Khairy M, Abbas W, Megahed AM, Emam MS. Semi analytical solution of MHD and heat transfer of couple stress fluid over a stretching sheet with radiation in porous medium. Front Heat Mass Transf. 2025;23(6):1833–46. doi:10.32604/fhmt.2025.069711. [Google Scholar] [CrossRef]

19. Ledinegg M. Instability of flow during natural and forced circulation. Die Wärme. 1938;61(48):891–8. [Google Scholar]

20. Bai S, Ma C, Peng X, Wen S. Thermoelastohydrodynamic behavior of gas spiral groove face seals operating at high pressure and speed. J Tribol. 2015;137(2):021502. doi:10.1115/1.4029358. [Google Scholar] [CrossRef]

21. Bergles AE, Kandlikar SG. On the nature of critical heat flux in microchannels. J Heat Transf. 2005;127(1):101–7. doi:10.1115/1.1839587. [Google Scholar] [CrossRef]

22. Kuo CJ, Peles Y. Pressure effects on flow boiling instabilities in parallel microchannels. Int J Heat Mass Transf. 2009;52(1–2):271–80. doi:10.1016/j.ijheatmasstransfer.2008.06.015. [Google Scholar] [CrossRef]

23. Odom BA, Miner MJ, Ortiz CA, Sherbeck JA, Prasher RS, Phelan PE. Microchannel two-phase flow oscillation control with an adjustable inlet orifice. J Heat Transf. 2012;134(12):122901. doi:10.1115/1.4007202. [Google Scholar] [CrossRef]

24. Barzegar Gerdroodbary M, Mousavi SV, Abdollahi SA. Influence of multiple electromagnetic sources for heat transfer improvement of ferrofluid flow inside the serpentine tube: a computational study. Comput Model Eng Sci. 2026;146(2):1–10. doi:10.32604/cmes.2026.076115. [Google Scholar] [CrossRef]

25. Ma X, Hu C, Chen H, Cu W, Zhang Y, Yang X, et al. Saturated/subcooled flow boiling heat transfer characteristics for R134a, HFE-7100 and Deionized water inside micro/mini-channels: part I—visualization of flow patterns and new diabatic two-phase flow pattern transition regimes. Int J Heat Mass Transf. 2024;234(7–8):126077. doi:10.1016/j.ijheatmasstransfer.2024.126077. [Google Scholar] [CrossRef]

26. Menter FR. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 1994;32(8):1598–605. doi:10.2514/3.12149. [Google Scholar] [CrossRef]

27. Da Riva E, Del Col D, Garimella SV, Cavallini A. The importance of turbulence during condensation in a horizontal circular minichannel. Int J Heat Mass Transf. 2012;55(13–14):3470–81. doi:10.1016/j.ijheatmasstransfer.2012.02.026. [Google Scholar] [CrossRef]

28. Zhang J, Li W, Minkowycz WJ. Numerical simulation of condensation for R410A at a different saturation temperature in mini/micro tubes. Numer Heat Transf Part A Appl. 2016;69(8):825–40. doi:10.1080/10407782.2015.1090844. [Google Scholar] [CrossRef]

29. Han J, Liu Y, Chu W, Zhao C, Bo H. Experimental study on visualized flow boiling in a narrow rectangular channel. Int Commun Heat Mass Transf. 2022;138(1):106383. doi:10.1016/j.icheatmasstransfer.2022.106383. [Google Scholar] [CrossRef]

30. Khovalyg D. Instability of flow boiling in parallel microchannels [dissertation]. Champaign, IL, USA: University of Illinois at Urbana-Champaign; 2017. [Google Scholar]