© The Author(s) 2025. Published by Science Press and EDP Sciences.

INTRODUCTION

The occurrence, development, and dissipation of bubbles within two-phase flow phenomena are prevalent in both natural environments and a variety of industrial applications, as illustrated in Figure 1. Examples include the processes of water boiling and condensation, chemical engineering practices, and nuclear systems [1–3]. Bubbles serve as versatile carriers with significant potential for applications such as enhancing boiling heat transfer [4], targeted cancer treatment [5], and controlling nanofluidic transport [6]. Furthermore, in fluid physics, the microscopic dynamics of multiscale bubbles play a crucial role in shaping the macroscopic characteristics of two-phase flow, such as interphase energy-mass transfer area [7], two-phase distributions in the flow field [8], bubble interactions, and transitions of flow patterns [9]. Consequently, the distribution of bubble sizes is of paramount importance. Previous studies have shown that bubbles of varying sizes exhibit distinct dynamic characteristics [7,8,10–13]. As an illustration, consider the two-phase flow in a small pipe with an inner diameter of 30 mm at room temperature (25°C) and normal pressure (0.1 MPa). Under low gas content, the gas phase is characterized by small, spherical bubbles that may exhibit distortion. As pressure decreases and gas content increases, the coalescence and expansion of these small bubbles facilitate the emergence of cap-shaped and slug bubbles throughout the channel’s cross-section. Furthermore, the lift force [8] influences the migration of small discrete bubbles toward the low-velocity region and the movement of large cap-shaped and slug bubbles toward the high-velocity region in a continuous phase flow field. The distribution characteristics of bubbles of varying sizes at the micro-scale indicate the macroscopic phase distribution characteristics under different flow patterns. Meanwhile, bubbles of different sizes exhibit varying rising velocities. Combining these factors, the interactions between multi-scale bubbles and their impact on the surrounding flow field vary, with large bubbles typically creating low-pressure wake regions behind them [14]. In these wake regions, smaller bubbles may accelerate and merge with larger bubbles, leading to a dynamic interfacial structure that contributes significantly to the complexity of two-phase flow evolution.

Figure 1 Bubble dynamic behavior with flow evolution and its potential occurrence scenarios.

Notably, in industrial systems involving two-phase flow, the mass, momentum, and energy transfer of the two phases occur through the phase-interface structure [15]. Accurate characterization of the phase-interface structure in such systems is crucial for optimizing system design and enhancing system security. As discussed previously, the bubble size essentially influences the magnitude of the phase-interface structure [14]. For instance, in a two-phase flow system with consistent gas phase content, the small bubbles’ presence results in a larger phase interface area within the unit spatial control volume (Interfacial Area Concentration, IAC). Conversely, larger bubbles lead to a smaller phase interface area. To predict the evolution of the phase interface, researchers developed the interfacial area transport equation (IATE) [16]. IATE takes IAC as the transport object, dynamically predicting the IAC evolution by considering the bubbles’ coalescence, break-up, evaporation, condensation, etc. [17]. Initially, Wu et al. proposed an IATE model for the bubbly flow. Subsequently, Fu [18] and Sun [19] extended the IATE model further to the slug flow and churn-turbulent flow by dividing the bubbles into two groups to establish bubble-interaction models (the group-1 bubbles G₁, referring to small and distorted bubbles; the group-2 bubbles G₂, referring to cap and slug bubbles). Up to now, IATE has been further extended from the simple pipe and rectangular to the complex rod bundle channel [20]. Recently, a growing concern has been that the current IATE model is unable to give satisfactory predictions in the transition from the bubbly flow to the slug flow [17]. To overcome such a flaw of the current IATE, researchers have developed wake entrainment models combined with a hypothetical bubble size distribution (BSD) [21,22]. Essentially, the IATE model flaw described above is because the model was developed based on the flat number density distribution of bubbles [14]. In the flow-pattern transition, BSD varies dramatically, which leads to drastic changes in the IAC and void fraction. Consequently, precise micro-scale bubble size distribution evolution (BSDE) prediction is crucial for accurately characterizing the IAC and void fraction evolution, which is also meaningful in describing other two-phase flow behaviors in both natural and industrial settings.

This study introduces a mathematical model, known as the BSDE model, which accounts for the transportation of bubble particles. The bubbles are categorized into two groups like that in the two-group IATE, to address the multiscale nature of bubbles in the BSDE model. The model development incorporates considerations of pressure variation, as well as bubble coalescence and break-up to capture their effects on the BSDE. To characterize bubble interactions, the bubble coalescence due to random collision (RC) and wake entrainment and the bubble break-up due to turbulent vortex and flow-field shearing-off are considered in the BSDE model. The model’s performance is assessed using BSD data obtained within a tight lattice bundle channel under gas-liquid flow conditions. The results show that the BSDE model can accurately predict the evolution of bubble sizes, with a deviation of 5% for one-group flow and 15% for two-group flow.

FRAMEWORK OF THE GENERAL BSDE MODEL

General remarks

For the two-phase flow in a general channel geometry, the bubble interactions affect the evolution of the bubble size. In addition, when considering the gravity head and friction resistance, the pressure gradually decreases along the axial direction, resulting in an increase in the bubble size along the channel. According to the upstream BSD, the downstream BSD can be predicted by considering the bubble coalescence/break-up and the pressure variation along the channel, which is the principle of the BSDE model. Given significant differences in the hydraulic behavior and interactions between the small and large bubbles, the bubbles were divided into two groups according to the critical size, as denoted by Equation (1), to consider the bubble interaction between the different scales. The feasibility of the two-group bubble division has been widely verified in the two-group IATE model [21–23]. As a theoretical derivation, the framework of the general BSDE model in this section is derived based on a vertical gas-liquid two-phase flow in a small- or medium-sized pipe under a 1-D steady state. However, it should be noted that the proposed BSDE model can be applied to other geometric flow channels theoretically whose hydraulic diameter is equivalent to that of small- or medium-sized pipe, after combining the geometry characteristics of the specific flow channel.

$\begin{array}{c}{D}_{\text{cl}}=4\sqrt{\frac{\sigma }{g∆\rho }}.\end{array}$(1)

Governing equation of the PDF for BSD

As mentioned above, the current BSDE model was developed by considering the effect of pressure and bubble coalescence/break-up during the flow evolution. The schematic of the 1-D control volume of a vertical channel is shown in Figure 2. The upstream face of the control volume is marked as the i-th face, and the downstream face is marked as the i+1-th face. The cross-sectional area of the channel is represented as A (m²); the axial length of the control volume is Δz (m); the void fraction at the i-th face is α(i); the velocity of the bubbles with diameter D at the i-th face is v_g(i, D) (m/s); the proportion of bubbles within diameter (D − 0.5ΔD, D + 0.5ΔD) at the i-th face is f(i, D) × ΔD. The symbol ΔD represents the differential of D. f(i, D) is the PDF of BSD at the i-th face, indicating that, in the unit space, the proportion of bubbles located in different diameter intervals in the volume of the total bubbles:

Figure 2 One-dimensional control volume.

$\begin{array}{c}f(i,D)=\frac{{\displaystyle {\int }_{D-0.5∆D}^{D+0.5∆D}\frac{\text{π}}{6}{D}^{3}N\left(i,D\right)\text{d}D}}{∆D\times {\displaystyle {\int }_{D_{\text{min}}}^{D_{\text{max}}}\frac{\text{π}}{6}{D}^{3}N\left(i,D\right)\text{d}D}}.\end{array}$(2)

During time Δt, the total volume of bubbles within (D − 0.5ΔD, D + 0.5ΔD) passing through the i-th face is marked as G_i:

$\begin{array}{c}{G}_i=f(i,D)\times ∆D\times {\alpha }_i\times v_{\text{g}}\left(i,D\right)\times A\times ∆t.\end{array}$(3)

Further, during Δt, the total volume of bubbles within (D − 0.5ΔD, D + 0.5ΔD) passing through the i+1th face is marked as G_i+1 (m³):

$\begin{array}{c}{G}_{i+1}=f(i+1,D)\times ∆D\times {\alpha }_{i+1}\times v_{\text{g}}\left(i+1,D\right)\times A\times ∆t.\end{array}$(4)

With pressure variation, the bubbles will expand or contract. From the i-th face to the i+1-th face, the increase in the total volume of bubbles within (D − 0.5ΔD, D + 0.5ΔD) due to the pressure effect is denoted as G_p (m³):

$\begin{array}{c}{G}_{\text{p}}=f\left(i,{\left(\frac{P_{i+1}}{P_i}\right)}^{1/3}D\right){\left(\frac{P_{i+1}}{P_i}\right)}^{1/3}∆\left(D\right){\alpha }_i{v}_{\text{g}}\left(i,{\left(\frac{P_{i+1}}{P_i}\right)}^{1/3}D\right)A∆t\frac{P_i}{P_{i+1}}-f\left(i,D\right)∆\left(D\right){\alpha }_i{v}_{\text{g}}\left(i,D\right)A∆t,\end{array}$(5)

where P_i and P_i+1 are the absolute pressures at the i-th and i+1-th faces, respectively. v_g is an independent parameter, similar to IATE model [24], which can be obtained by experiment or from the momentum if the current model is combined with the two-fluid model. In addition, the introduction of the velocity indicates that Equation (5) also contains the velocity effect on BSD evolution.

With bubble coalescence and break-up, the increase in the total volume of bubbles within (D − 0.5ΔD, D + 0.5ΔD) is denoted as G_S (m³):

$\begin{array}{c}{G}_{\text{S}}=∆z\times A\times ∆t\times {\alpha }_i\times ∆\left(D\right){\displaystyle \sum _j^k{S}_j^i}\left(D\right),\end{array}$(6)

where $S_j^i\left(D\right)$ is the BSD source term for the j-th bubble interaction mechanism at the i-th face (1/(m s)).

The conservation equation of bubble volume concerning bubble size of (D − 0.5ΔD, D + 0.5ΔD) in the control volume can be written as

$\begin{array}{c}{G}_{i+1}-G_i=G_{\text{P}}+G_{\text{S}}.\end{array}$(7)

Rearrange Equation (7) as

$\begin{array}{c}f\left(i+1,D\right)=f\left(i,{\left(\frac{P_{i+1}}{P_i}\right)}^{1/3}\times D\right)\frac{{\alpha }_i}{{\alpha }_{i+1}}\frac{v_{\text{g}}\left(i,{\left(\frac{P_{i+1}}{P_i}\right)}^{1/3}D\right)}{v_{\text{g}}\left(i+1,D\right)}{\left(\frac{P_i}{P_{i+1}}\right)}^{2/3}\text{+}\frac{∆z}{v_{\text{g}}\left(i,D\right)}\times {\displaystyle \sum _j^k{S}_j^i}\left(D\right).\end{array}$(8)

In Equation (8), the LHS (Left-hand Side) is the BSD at the downstream node, and the RHS (Right-hand Side) represents the function of the upstream BSD. The first term on the RHS is the BSD source for the variation in pressure and bubble velocity, and the second term denotes the BSD source for the bubble coalescence and break-up. Equation (8) clearly shows the principle of the BSDE model, based on the BSD information at the upstream node, considering the effects of bubble velocity change, pressure change, and bubble interactions on BSD, to predict the BSD downstream.

Source terms for bubble interactions

For a realistic two-phase flow, the mechanisms of bubble coalescence and breakup can be divided into five categories, which affect the BSD evolutions, as shown in Figure 3: (a) bubble coalescence due to RC driven by turbulent vortices; (b) absorption of the trailing bubble by the leading bubble owing to wake entrainment (WE); (c) bubble breakup due to turbulent impact (TI); (d) shearing off (SO) of a large number of small bubbles from the tail of a large bubble; (e) large bubble breakup due to surface instability (SI). Given the different effects of channel size (hydraulic diameter) on multi-scale bubbles, the above bubble interaction mechanism does not usually occur simultaneously. For the two-phase flow in the small- and medium-sized pipes concerned in this study, unimportant bubble interaction mechanisms should be excluded during the development of the model. Specifically, in the small- and medium-sized channels, the G₂ bubbles’ size is in the same order of magnitude as the channel cross-section. The random movement of G₂ bubbles driven by the turbulent vortex is not significant. In addition, the RC sub-mechanisms that involved G₂ bubbles can be ignored (RC^{(12, 2)} and RC⁽²⁾). According to experimental observations, in the medium-size channel, the probability of G₂ bubbles breaking up owing to TI and SI is small; in addition, the breakup mechanisms of G₂ bubbles can be ignored (TI^{(2, 11)}, TI^{(2, 12)}, TI⁽²⁾, SI⁽²⁾). The remaining mechanisms were simulated: RC⁽¹⁾, RC^{(11, 2)}, WE⁽¹⁾, WE^{(11, 2)}, WE^{(12, 2)}, WE⁽²⁾, TI⁽¹⁾, and SO^{(2, 12)} (Supplementary information). However, given the complexity of the real two-phase flow, the neglected source/sink terms may also have important effects on the BSD evolution under different flow conditions and channel geometry. In addition, the current development of the source/sink model assumes the range of bubble sizes applicable to different models, which is also a source of error. Nevertheless, the current data evaluation results prove the rationality and feasibility of the current modeling strategy. In future studies, more data from different flow conditions and channel geometries will be needed to verify and optimize the model’s performance.

Figure 3 Mechanisms of bubble coalescence and breakup. (a) RC; (b) WE; (c) TI; (d) SO; (e) SI.

Derivation method of undetermined coefficients and BSDE model assessment criteria

During the development of the BSDE model, undetermined coefficients were introduced in the submodel of each source term, including C_{RC1_1} in the G₁ bubble coalescence due to RC, C_{WE1_1} in the G₁ bubble coalescence due to wake entrainment, C_{WE2_1} in the G₂ bubble coalescence due to wake entrainment, C_{WE21_1} in the G₁ and G₂ bubble coalescence due to wake entrainment, C_{SO21_1} in the G₂ bubble SO, and C_{TI1_1} in the G₁ bubble breakup due to TI. After the undetermined coefficients are determined based on the bubble interaction dominant experimental results, the BSDE model can be experimentally evaluated. The schematic of the computational scheme of the BSDE model is shown in Figure 4.

Figure 4 Schematic of the computational scheme of the BSDE model.

To quantify the differences between the experimental results and the theoretical predictions, the relative error was determined by the area surrounded by the experimentally measured BSD, the model-predicted BSD, and the x-axis, as shown in Figure 5. In Figure 5, a larger shadowed area would imply a more accurate model. The expression of the model relative error is

Figure 5 Relative error in BSD prediction.

$\begin{array}{c}\text{errorBSD}=\frac{1}{2}{\displaystyle {\sum }_k\left[1-{\displaystyle {\int }_{D_{\min }}^{D_{\max }}\min \left(f_{k,\text{pred}}\left(D\right),f_{k,\exp }\left(D\right)\right)\text{d}D}\right]},\end{array}$(9)

where f_k,pred(D) is the BSD from the model prediction at the k-th port, and f_k,exp(D) is the experimentally obtained BSD at the k-th port.

APPLICATION OF THE BSDE MODEL TO A TIGHT-LATTICE ROD BUNDLE AND EVALUATION

In this section, the BSDE model proposed in Section “Framework of the general BSDE model” was applied to predict the BSD evolution in a specific tight-lattice bundle. The experimental database, model coefficients derivation, and prediction results were discussed in the following.

Experimental database and coefficients derivation

The experimental data were from a dual sub-channel tight lattice rod bundle, whose flow conditions are shown in Figure 6, including the bubbly flow, cap-bubbly flow, and slug flow. The experiments used a dual wire-mesh sensor to measure the phase interface information with a high resolution shown in Figure 7 and to acquire the void fraction, BSD, IAC, and bubble velocity with the Maxwell No Cut method [25] and self-developed Eulerian post-processing algorithm [26]. Figure 8 shows the typical BSD evolution experimental results. The pressure data at the measurement ports were also recorded through differential pressure transmitters. More details can be obtained from the previous research [13]. By linear interpolation, the flow parameters, including the group-wise bubble velocity, void fraction, and absolute pressure along the channel axial direction, were calculated based on the experimental data. This avoids the effects of model uncertainty in predicting these parameters on the BSDE model. The BSDs measured in the experiments at z/D_h = 57.90 were used as the inlet boundary of the BSDE model, and the experimental data at z/D_h = 86.86 and 115.8 were used for fitting the BSDE model coefficients and evaluating the model. The calculation method of BSD can be referred to from Refs. [26,27].

Figure 6 Experimental flow conditions at z/Dh = 115.8.

Figure 7 3-D phase interface distribution in the tight lattice bundle based on the WMS technique. (a) Bubbly flow; (b) cap-bubbly flow; (c) slug flow [28].

Figure 8 Typical BSD evolution experimental results. (a) One-group flow, case 20; (b) two-group flow, case 11.

The model coefficient is determined according to the dominant flow condition of the mechanism. It is worth mentioning that the BSDE model requires the geometrical information of G₂ bubbles to calculate the BSDE source term from the SO of a G₂ bubble (SO⁽²¹⁾). Given the absence of modeling of bubble geometric characteristics in the current study, to ensure the integrity of the model, the fitting correlations obtained from the experimental data were used [29] to calculate the G-2 bubbles’ height and projected area diameter, as shown in Equations (10) and (11). This does not affect the overall model framework’s valuation and the source/sink terms’ functionality evaluation. These experimental data were also used for the BSDE-model evaluation in Sections “Evaluation in one-group flow conditions” and “Evaluation in two-group flow conditions”.

$\begin{array}{c}H\left(D_{\text{3}}\right)=-\text{0.004+20.3}\times D_{\text{3}}{}^{\text{1.566}},\end{array}$(10)

$\begin{array}{c}{D}_{\text{bd}}\left(D_3\right)=-\text{0.00568+0.11}\times D_{\text{3}}{}^{\text{0.426}}.\end{array}$(11)

The following is the coefficient fitting process. It follows the IATE model coefficient confirmation strategy, where the experimental data for fitting the undetermined coefficients and evaluating the model were not distinguished [21,30,31]. First, the experimental data for one-group flow were used to fit the coefficients in the G₁ bubble interaction mechanisms. Six flowing conditions of one-group flow, with different superficial liquid velocities, were selected to fit these coefficients. In cases 21, 32, and 35, the bubble size increases along the channel owing to bubble coalescence. These three conditions were selected to fit C_{RC1_1}, with coefficients C_{TI1_1} and C_{WE1_1} set to zero. In cases 8, 16, and 30, the effects of bubble coalescence on global BSD were compensated by the bubble breakup due to TI. These three flowing conditions were used for C_{TI1_1}, with C_{WE1_1} set to zero. Although the source term distribution of WE⁽¹⁾ was developed, WE⁽¹⁾ has a negligible influence on the experimental data of BSDE. Therefore, the effects of WE⁽¹⁾ were incorporated into RC⁽¹⁾, and the coefficient C_{WE1_1} was set to zero. After C_{TI1_1} was fitted, C_{RC1_1} was refitted with the updated C_{TI1_1} and into the next fitting loop until the values of C_{RC1_1} and C_{TI1_1} were stable.

Second, the coefficients of G₂ bubble interactions were fitted after C_{RC1_1} and C_{TI1_1} were determined. Cases 1, 12, and 23 were selected to fit C_{WE2_1}, with C_{SO21_1} and C_{WE21_1} set to zero. In these three conditions, the size of G₂ bubbles was significantly larger owing to WE⁽²⁾. In cases 11 and 24, a large number of G₁ bubbles were generated by SO⁽²¹⁾; therefore, these two flowing conditions are suitable for C_{SO21_1} fitting. Coefficient C_{WE21_1} was set to zero during C_{SO21_1} fitting. In cases 1 and 12, the generation of G₁ bubbles from SO⁽²¹⁾ was compensated by the consumption of G₁ bubbles from WE⁽²¹⁾. These two flowing conditions were used for C_{WE21_1} fitting. After C_{WE21_1} was fitted, C_{WE2_1} was refitted with the updated C_{SO21_1} and C_{WE21_1} and into the next fitting loop until the values of these three coefficients were stable.

The fitted coefficients are shown in Table 1.

Evaluation in one-group flow conditions

The relative error of the BSDE model for one-group flow (bubbly flow) is shown in Figure 9. All the flowing conditions display a relative error smaller than 5%, indicating that the BSDE model has high accuracy for BSD prediction in the one-group flow. It is worth noting that the evolution of BSD in the upstream and downstream is insignificant, which is another reason for such a small deviation. The detailed BSD prediction results for cases 33, 18, and 6 are shown in Figures 10–12, respectively. Figures 10(a), 11(a), and 12(a) show the experimentally measured BSDs at different axial positions of the flow channel. Figure 10(b) and (c), Figure 11(b) and (c), and Figure 12(b) and (c) show the BSD comparison between the predicted and experimental measurements. Moreover, Figures 10(d), 11(d), and 12(d) show the BSD accumulation change caused by each source/sink term under different bubble diameters during the whole flow development process.

Figure 9 Model evaluation of one-group flow.

Figure 10 BSDE model evaluation under case 33, with jf = 1.86 m/s, jg = 0.11 m/s. (a) Experimental BSD; (b) model prediction vs. experiments for z/Dh = 86.86; (c) model prediction vs. experiments for z/Dh = 115.8; (d) accumulation of each source term along the channel length.

Figure 11 BSDE model evaluation under case 18, with jf = 1.39 m/s and jg = 0.16 m/s. (a) Experimental BSD; (b) model prediction vs. experiments for z/Dh = 86.86; (c) model prediction vs. experiments for z/Dh = 115.8; (d) accumulation of each source term along the channel length.

Figure 12 BSDE model evaluation under case 6, with jf = 0.93 m/s and jg = 0.17 m/s. (a) Experimental BSD; (b) model prediction vs. experiments for z/Dh = 86.86; (c) model prediction vs. experiments for z/Dh = 115.8; (d) accumulation of each source term along the channel length.

For case 33, the peak value of the BSD at the first measured port (z/D_h = 57.90) appears at D_b = 2.5 mm. With the flow development, the bubble size increases and the BSD profile shifts toward the right side. The peak value of the BSD at the third port (z/D_h = 115.8) appears at 2.8 mm, as shown in Figure 10(a). The comparison of the BSDs between the model prediction and experimental measurements for z/D_h = 86.86 and 115.8 are shown in Figure 10(b) and (c), respectively, displaying good consistency.

In the one-group flow, the bubble expansion and RC⁽¹⁾ are the dominant mechanisms for the BSD evolution. The bubble expansion shifts the BSD toward the right side, which in the source term distribution consumes small-sized bubbles and generates large-sized bubbles, which also implies the increment of the IAC. Figure 6(d) shows that bubbles with D_b = ~2 mm were consumed and those with D_b = ~4 mm were generated in RC⁽¹⁾. In addition, TI⁽¹⁾ consumed large bubbles and generated small bubbles. In case 33, with j_f = 1.86 m/s, the fluctuation velocity is large enough for the breakup of bubbles with D_b = 4–5 mm, as shown in Figure 10(d), and generates bubbles with D_b = 3–4 mm. With the decrease in the superficial liquid velocity, the velocity of fluctuation decreases, the diameter of the broken bubbles increases, and the contribution of TI⁽¹⁾ decreases. For case 18, with j_f = 1.39 m/s, only a few bubbles with diameters of approximately 5–6 mm are consumed due to TI⁽¹⁾ and that of generated bubbles is D_b = 4–5 mm. For case 6, with j_f = 0.93 m/s, the contribution of TI⁽¹⁾ can be ignored. As the BSD profile does not display perfect smoothness, the source term accumulation and distribution also show a certain degree of fluctuation.

Evaluation in two-group flow conditions

The relative error of the BSDE model for the prediction of a two-group flow is shown in Figure 13. All the flowing conditions show a relative error smaller than 15%, indicating that the BSDE model accurately predicts the BSDE along the channel. In a two-group flow, the effects of the G₂ bubble coalescence due to WE⁽²⁾ and SO on BSDE were significant. The BSDE was more significant for a two-group flow than for a one-group flow. The detailed BSD prediction results for cases 23, 10, and 1 are shown in Figures 14–16, respectively. Figures 14(a), 15(a), and 16(a) show the experimentally measured BSD at different axial positions of the flow channel. Figure 14(b) and (c), Figure 15(b) and (c), and Figure 16(b) and (c) show the BSD comparison between the model prediction and experimental measurements, respectively. Further, similar to that in the one-group flow conditions, Figure 14(d) and (e), Figure 15(d) and (e), and Figure 16(d) and (e) show separately the G₁ and G₂ BSD accumulation change caused by each source/sink term under different bubble diameters during the whole flow development process.

Figure 13 Model evaluation of a two-group flow.

Figure 14 BSDE model evaluation under case 23, with jf = 1.86 m/s, jg = 0.70 m/s. (a) Experimental BSD; (b) model prediction vs. experiments for z/Dh = 86.86; (c) model prediction vs. experiments for z/Dh = 115.8; (d) accumulation of each source term on G1 bubbles along the channel length; (e) accumulation of each source term on G2 bubbles along the channel length.

Figure 15 BSDE model evaluation under case 10, with jf = 1.39 m/s, jg = 0.68 m/s. (a) Experimental BSD; (b) model prediction vs. experiments for z/Dh = 86.86; (c) model prediction vs. experiments for z/Dh = 115.8; (d) accumulation of each source term on G1 bubbles along the channel length; (e) accumulation of each source term on G2 bubbles along the channel length.

Figure 16 BSDE model evaluation under case 1, with jf = 0.93 m/s, jg = 0.49 m/s. (a) Experimental BSD; (b) model prediction vs. experiments for z/Dh = 86.86; (c) model prediction vs. experiments for z/Dh = 115.8; (d) accumulation of each source term on G1 bubbles along the channel length; (e) accumulation of each source term on G2 bubbles along the channel length.

In a two-group flow, the bubbles were naturally divided into a small-size group and a large-size group, as shown in Figure 14(a). For example, consider case 23, in which the G₁ bubbles show D_b < 11 mm and the G₂ bubbles show 20 < D_b < 40 mm. The range of 11 < D_b < 20 mm is the transition region with a few bubbles. With the flow development, the BSD profile of the G₂ bubbles shifts toward the right side owing to WE⁽²⁾, and the BSD profile of the G₁ bubbles increases owing to SO⁽²¹⁾. Therefore, the G₂ IAC will decrease, and the G₁ IAC will increase. In addition, the profile of G₁ bubble BSD remains similar at different axial positions, proving the rationality of the assumption of the generated BSD concerning SO⁽²¹⁾, as discussed in Supplementary Section 2.5. It was assumed that the BSD of the generated G₁ bubbles from SO is proportional to the BSD of the G₁ bubble at the current node.

The comparisons of BSD between the model prediction and experimental measurements for z/D_h = 86.86 and 115.8 are shown in Figure 14(b) and (c), respectively, displaying good consistency. The BSD variation characteristics of the bubbles in the two groups are reproduced by the BSDE model. The accumulation of each source term on the G₁ bubbles is shown in Figure 14(d). Multiple bubble interaction mechanisms were activated and the SO⁽²¹⁾, TI1, Expa, and RC1 are observed as the dominant mechanism for the G₁-bubble BSDE. SO⁽²¹⁾ is an intergroup mechanism, which consumed G₂ bubbles and generated G₁ bubbles. Although SO⁽²¹⁾ is a dominant source term for the G₁ bubble, it has a negligible effect on the G₂ bubbles, as shown in Figure 14(e). For the G₂ bubbles, the WE⁽²⁾ and Expa is the dominant mechanism, which is consistent with the experimental observations. In WE⁽²⁾, bubbles with 25 mm < D_b < 35 mm were consumed, and bubbles with 35 mm < D_b < 45 mm were generated.

In case 10, with j_f = 1.39 m/s and j_g = 0.68 m/s, a significant peak was observed in the BSD of the G₂ bubble for z/D_h = 57.90 (as shown by the green line in Figure 15(a)). Most of the G₂ bubbles had diameters ranging from 28 to 38 mm. These were observed under the flowing conditions near the flow regime transition boundary. Many G₁ bubbles were generated at the inlet, and then drastic bubble coalescence occurred, resulting in small G₂ bubbles with similar sizes. With the flow development, G₂ bubbles coalesce owing to WE⁽²⁾, and the BSD of the G₂ bubble flattens further. The BSDE of the G₂ bubbles was accurately predicted by the BSDE model, as shown in Figure 15(b) and (c). WE⁽²⁾ was also the dominant mechanism for the G₂-bubble BSDE, as shown in Figure 15(e). Bubbles with a diameter of 30–35 mm were consumed and bubbles with a diameter of 35–50 mm were generated. For the G₁-bubble BSDE, the SO⁽²¹⁾ was the dominant mechanism. The comparison results for case 1 are similar to those of cases 23 and 10, as shown in Figure 16. Thus, the BSDE characteristics were accurately predicted by the BSDE model. The SO⁽²¹⁾ was the dominant mechanism for the G₁-bubble BSDE and WE⁽²⁾ was the dominant mechanism for the G₂ bubbles. With the decrease in superficial liquid velocity, the fluctuation amount of the liquid velocity decreases, and the effects of TI⁽¹⁾ on the BSDE decrease. Therefore, the contribution of TI⁽¹⁾ was significant in case 23 with j_f = 1.86 m/s, whereas it was weak in case 1 with j_f = 0.93 m/s.

CONCLUSION

The BSD is a critical parameter in gas-liquid two-phase flow, significantly influencing characteristics such as void fraction and interfacial area. Pressure fluctuations and bubble interactions during flow progression result in BSDE along the channel. This study introduces a 1-D BSDE model to simulate the evolution of BSD along a vertical channel.

Due to the complex multiscale nature of the two-phase flow, bubbles are categorized into two groups to model the distribution of the bubble interaction terms. Six distinct bubble interaction mechanisms were developed in this study: (1) coalescence of G₁ bubbles due to RC; (2) coalescence of G₁ bubbles due to wake entrainment; (3) coalescence of G₂ bubbles due to wake entrainment; (4) coalescence of G₂ and G₁ bubbles due to wake entrainment; (5) SO of G₂ bubbles; and (6) G₁ bubble breakup due to TI. Each mechanism encompasses the distributions of bubble consumption and production. Through coefficient constraints, it was determined that the total volume of consumed bubbles equaled the total volume of generated bubbles, ensuring conservation of bubble volume.

Experimental data from an air-water two-phase flow in a double sub-channel tight lattice rod bundle were utilized to calibrate the BSDE model and assess its precision. The flow region includes one-group bubbly flow and two-group cap-bubbly flow. The model inlet boundary was set at z/D_h = 57.90, with validation data obtained at z/D_h = 86.86 and 115.8. In the case of one-group flow, the BSD exhibited a single-peak profile with bubble sizes smaller than 11 mm, and the BSDE model’s prediction deviation was within 5%. The dominant mechanisms influencing the BSD evolution were identified as bubble expansion and RC⁽¹⁾, leading to a shift of the BSD towards higher values. In the case of two-group flow, the BSDE model’s prediction deviation was within 15%. The effects of the G₂ bubble coalescence due to WE⁽²⁾ and the G₂ bubble breakup due to SO on BSDE were significant. With the flow development, the BSD profile of the G₂ bubbles shifts toward the right side owing to WE⁽²⁾, and the BSD profile of G₁ bubbles increases owing to SO⁽²¹⁾. The accurate predictions of BSDE demonstrate the feasibility of the proposed BSDE model.

The verification of the model presented in this study is based on data from a tight lattice rod bundle. It is worth noting that, the source/sink models of the BSDE model were developed based on the small- and medium-sized pipes. It implies that the current model can theoretically be extended to geometric channels of different cross-section shapes, as long as their hydraulic diameters are in the small- and medium-sized pipe category, with minor adjustments, as only the G₂ bubble-shape assumption and SO mechanism are related to channel shape. Future research should evaluate the BSDE model’s accuracy using experimental data in basic geometric channels like circular or rectangular ones. Additionally, the void fraction distribution prediction model and IATE model should be improved based on the BSDE model in future studies.

Funding

This work was supported by the National Natural Science Foundation of China (12322510 and 12275174), the Shanghai Rising-Star Program (22QA1404500), the Science and Technology Commission of Shanghai Municipality (24DZ3100300), and the Lingchuang Project of China National Nuclear Corporation.

Author contributions

In this project, Y.X. was responsible for conceptualization, methodology, investigation, resources, writing review & editing, supervision, and funding acquisition; H.Z. was responsible for conceptualization, software, methodology, investigation, data curation, and writing original draft; X.Y. was responsible for writing review & editing, writing original draft, software, and investigation; H.G. was responsible for investigation, validation, resources, and funding acquisition.

Conflict of interest

The authors declare no conflict of interest.

Supplementary information

Supplementary file provided by the authors. Access here

The supporting materials are published as submitted, without typesetting or editing. The responsibility for scientific accuracy and content remains entirely with the authors.

References

All Tables

All Figures

	Figure 1 Bubble dynamic behavior with flow evolution and its potential occurrence scenarios.
In the text

	Figure 2 One-dimensional control volume.
In the text

	Figure 3 Mechanisms of bubble coalescence and breakup. (a) RC; (b) WE; (c) TI; (d) SO; (e) SI.
In the text

	Figure 4 Schematic of the computational scheme of the BSDE model.
In the text

	Figure 5 Relative error in BSD prediction.
In the text

	Figure 6 Experimental flow conditions at z/Dh = 115.8.
In the text

	Figure 7 3-D phase interface distribution in the tight lattice bundle based on the WMS technique. (a) Bubbly flow; (b) cap-bubbly flow; (c) slug flow [28].
In the text

	Figure 8 Typical BSD evolution experimental results. (a) One-group flow, case 20; (b) two-group flow, case 11.
In the text

	Figure 9 Model evaluation of one-group flow.
In the text

	Figure 10 BSDE model evaluation under case 33, with jf = 1.86 m/s, jg = 0.11 m/s. (a) Experimental BSD; (b) model prediction vs. experiments for z/Dh = 86.86; (c) model prediction vs. experiments for z/Dh = 115.8; (d) accumulation of each source term along the channel length.
In the text

	Figure 11 BSDE model evaluation under case 18, with jf = 1.39 m/s and jg = 0.16 m/s. (a) Experimental BSD; (b) model prediction vs. experiments for z/Dh = 86.86; (c) model prediction vs. experiments for z/Dh = 115.8; (d) accumulation of each source term along the channel length.
In the text

	Figure 12 BSDE model evaluation under case 6, with jf = 0.93 m/s and jg = 0.17 m/s. (a) Experimental BSD; (b) model prediction vs. experiments for z/Dh = 86.86; (c) model prediction vs. experiments for z/Dh = 115.8; (d) accumulation of each source term along the channel length.
In the text

	Figure 13 Model evaluation of a two-group flow.
In the text

	Figure 14 BSDE model evaluation under case 23, with jf = 1.86 m/s, jg = 0.70 m/s. (a) Experimental BSD; (b) model prediction vs. experiments for z/Dh = 86.86; (c) model prediction vs. experiments for z/Dh = 115.8; (d) accumulation of each source term on G1 bubbles along the channel length; (e) accumulation of each source term on G2 bubbles along the channel length.
In the text

	Figure 15 BSDE model evaluation under case 10, with jf = 1.39 m/s, jg = 0.68 m/s. (a) Experimental BSD; (b) model prediction vs. experiments for z/Dh = 86.86; (c) model prediction vs. experiments for z/Dh = 115.8; (d) accumulation of each source term on G1 bubbles along the channel length; (e) accumulation of each source term on G2 bubbles along the channel length.
In the text

	Figure 16 BSDE model evaluation under case 1, with jf = 0.93 m/s, jg = 0.49 m/s. (a) Experimental BSD; (b) model prediction vs. experiments for z/Dh = 86.86; (c) model prediction vs. experiments for z/Dh = 115.8; (d) accumulation of each source term on G1 bubbles along the channel length; (e) accumulation of each source term on G2 bubbles along the channel length.
In the text