Experimental and Numerical Study of the Flow and Heat Transfer in a Bubbly Turbulent Flow in a Pipe with Sudden Expansion

The flow patterns and heat transfer of a downstream bubbly flow in a sudden pipe expansion are experimentally and numerically studied. Measurements of the bubble size were performed using shadow photography. Fluid phase velocities were measured using a PIV system. The numerical model was employed the Eulerian approach. The set of RANS equations was used for modelling two-phase bubbly flows. The turbulence of the carrier liquid phase was predicted using the Reynolds stress model. The peak of axial and radial fluctuations of the carrier fluid (liquid) velocity in the bubbly flow is observed in the shear layer. The addition of air bubbles resulted in a significant increase in the heat transfer rate (up to 300%). The main enhancement in heat transfer is observed after the point of flow reattachment.


Introduction
Gas-liquid bubbly flows are frequently used in many engineering and practical applications. These flows are usually turbulent with a strong interfacial coupling between the fluid carrier and dispersed phases. A study of such flows is usually complicated by flow separation at sharp edges, and heat transfer. Flow recirculation largely determines the mean and turbulent characteristics and has a significant influence on the transport processes [1,2]. The modeling of bubble motion and their distributions over a duct or pipe cross-section is very important for safety and for the simulation of various emergency scenarios in steam and nuclear power equipment.
Experimental works [3,4] were probably among the first studies of bubbly flows downstream of a sudden expansion. They studied bubbly flows behind a sudden expansion in a horizontal duct. The pressure drop, void fraction, bubble size, and the mean and fluctuation velocities of the phases, were measured in these papers. Aloui and Souhar [3] and Rinne and Loth [4] experimentally observed a reduction in diameter of bubbles along the duct length due to the growth in the longitudinal pressure gradient and the bubbles break up in separation zone.
Later experimental [5][6][7][8] and numerical [9,10] studies of separated isothermal bubbly flows were continued. A theoretical and experimental study of upward bubbly flow in a pipe with sudden expansion [5] was carried out on. The mean and fluctuating bubble velocities, local void fraction, bubble distributions, pressure drop and friction on the wall were measured. The distributions of the bubble's velocity fluctuations were similar to that of the carrier fluid phase. Measurements were performed in the horizontal flow of the air-oil mixture [6], which was transient in the annular flow remove excess heat from the tank, a heat exchanger, through which cooling water was pumped, was installed. The liquid temperature in the tank was measured by a resistive temperature sensor installed in the tank and automatically controlled by management of the cooling water flow rate by a type TRM-2M programmable controller. This makes it possible to keep the liquid temperature at the beginning of the test section in the 24.9−25.1 °C range, according to an analysis of the signal from the temperature sensor mounted in the test section's inlet, which is a thin film platinum sensor. The nominal resistance of this sensor is 1000 Ohm at 0 °C. The response time of this type of sensor in water flow is about 0.2 s. The measurement uncertainty of temperature measurement in the test section inlet is about ±0.15 °C. Figure 1. The scheme of the experimental facility. 1-main tank; 2-pump; 3-supply line; 4-control valve; 5-ultrasound flow meter; 6-gas-liquid mixer; 7-gas flow controller; 8-the pipe of small diameter (before expansion); 9-expansion area; 10-box filled by water; 11-video recorder; 12-test section; 13-copper conductors; 14-infrared (IR) camera; 15-resistive temperature detectors (RTDs); 16-the pipe of large diameter (after expansion); 17-downcomer. A Grundfos CHI 4-40 (Ql ≤ 4.5 m 3 /h) is used to pump the liquid. The liquid flux is directed upstream, the liquid flow rate is controlled using a valve, 4. The flow rate is monitored using an ultrasonic flow meter, 5. In this paper the effect of intermediate sized bubbles (d = 1 ÷ 3 mm) on the heat transfer in the pipe with sudden expansion is studied. Bubbles are generated by feeding air into the liquid flow through 12 capillaries, where the length of each capillary is 50 mm and inner diameter is 0.7 mm. The capillaries are uniformly distributed across the pipe channel's cross-section. The gas flow rate is determined using a gas flow controller, 7, produced by Bronkhorst (Qg ≤ 1 L/min). The measurement error of the gas and liquid flow rates is ±2% of the measured value. Then, a gas-liquid mixture is fed into the test section. It is a small-diameter pipe, 8 (2R1 = 15 mm), which is inserted through the adapter and into a large diameter pipe. The inner diameter of the pipe, 16, is 2R2 = 42 mm, and the step height is H = 13.5 mm. The expansion ratio is defined as ER = (R2/R1) 2 = 7.84. The end of the small pipe is mounted flush with the plane of the adapter located in the flow separation zone. To ensure that a fully developed gas-liquid flow is obtained in the inlet of the test section, the length of the flow-stabilization area before the measurement region is 140R1. The experimental setup and the method for heat transfer measurements was presented in [15] in details.

PIV/PLIF Measurements
In order to perform the fluid phase velocity distribution measurements, a particle image velocimetry (PIV) "Polis-PIV" system was used. The Plexiglass pipe was mounted instead of the heat transfer measurement unit for PIV measurements and it was installed into a rectangular box that was filled with water to reduce image distortions. In two-phase flow, bubbles illuminated by a laser light produce a bright glare that can damage the image sensor. That's why the planar laser-induced Figure 1. The scheme of the experimental facility. 1-main tank; 2-pump; 3-supply line; 4-control valve; 5-ultrasound flow meter; 6-gas-liquid mixer; 7-gas flow controller; 8-the pipe of small diameter (before expansion); 9-expansion area; 10-box filled by water; 11-video recorder; 12-test section; 13-copper conductors; 14-infrared (IR) camera; 15-resistive temperature detectors (RTDs); 16-the pipe of large diameter (after expansion); 17-downcomer.
A Grundfos CHI 4-40 (Ql ≤ 4.5 m 3 /h) is used to pump the liquid. The liquid flux is directed upstream, the liquid flow rate is controlled using a valve, 4. The flow rate is monitored using an ultrasonic flow meter, 5. In this paper the effect of intermediate sized bubbles (d = 1 ÷ 3 mm) on the heat transfer in the pipe with sudden expansion is studied. Bubbles are generated by feeding air into the liquid flow through 12 capillaries, where the length of each capillary is 50 mm and inner diameter is 0.7 mm. The capillaries are uniformly distributed across the pipe channel's cross-section. The gas flow rate is determined using a gas flow controller, 7, produced by Bronkhorst (Qg ≤ 1 L/min). The measurement error of the gas and liquid flow rates is ±2% of the measured value. Then, a gas-liquid mixture is fed into the test section. It is a small-diameter pipe, 8 (2R 1 = 15 mm), which is inserted through the adapter and into a large diameter pipe. The inner diameter of the pipe, 16, is 2R 2 = 42 mm, and the step height is H = 13.5 mm. The expansion ratio is defined as ER = (R 2 /R 1 ) 2 = 7.84. The end of the small pipe is mounted flush with the plane of the adapter located in the flow separation zone. To ensure that a fully developed gas-liquid flow is obtained in the inlet of the test section, the length of the flow-stabilization area before the measurement region is 140R 1 . The experimental setup and the method for heat transfer measurements was presented in [15] in details.

PIV/PLIF Measurements
In order to perform the fluid phase velocity distribution measurements, a particle image velocimetry (PIV) "Polis-PIV" system was used. The Plexiglass pipe was mounted instead of the heat transfer measurement unit for PIV measurements and it was installed into a rectangular box that was filled with water to reduce image distortions. In two-phase flow, bubbles illuminated by a laser light produce a bright glare that can damage the image sensor. That's why the planar laser-induced fluorescence (PLIF) approach was applied in our measurements [16]. Fluorescent seeding particles made of polymethyl methacrylate filled with Rhodamine B produced by Dantec Dynamics\(hydrophobic, size distribution 1-20 µm, wavelength range 550-700 nm) were used in our experiments. A threshold optical filter was used to prevent glare from bubbles surfaces on the final images. Thus, we obtained positions of tracers on the initial experimental data, the positions of the bubbles are invisible. This allows us to obtain the velocity field for liquid while ignoring the characteristics of the gas bubbles. This relates to complex shapes of bubbles, which are quite different from spherical or elliptical in the shear region of the flow (see Figure 2). The recognition procedure for such bubbles is a very complex problem.
Energies 2018, 11, x FOR PEER REVIEW 4 of 18 made of polymethyl methacrylate filled with Rhodamine B produced by Dantec Dynamics\(hydrophobic, size distribution 1-20 μm, wavelength range 550-700 nm) were used in our experiments. A threshold optical filter was used to prevent glare from bubbles surfaces on the final images. Thus, we obtained positions of tracers on the initial experimental data, the positions of the bubbles are invisible. This allows us to obtain the velocity field for liquid while ignoring the characteristics of the gas bubbles. This relates to complex shapes of bubbles, which are quite different from spherical or elliptical in the shear region of the flow (see Figure 2). The recognition procedure for such bubbles is a very complex problem. One thousand pairs of PIV images were obtained for both of two cases. One of these cases was a single-phase flow with ReH = 2.09 × 10 4 . In another case, a gas phase with β = 3.5% was added to the flow. Image processing was carried out using the "ActualFlow" software. At first a procedure of subtraction of the mean intensity field averaged over the whole sample range was applied for the raw data. The image processing was carried out by means of iterative cross-correlation algorithm. The overlap of interrogation windows was 75%. The size of the interrogation window was 64 × 64 pixels. Following validation procedures were used: peak validation with the threshold of 1.2, adaptive median 3 × 3 filter and range validation taking into account of maximum velocity on the axis of the small pipe multified by a coefficient 1.2. The uncertainty of the liquid velocity measurements was about 5%.

Bubble Size Measurements
Measurements of the sizes of bubbles are performed using video with shadow illumination prior to the location of the expansion. To reduce the optical distortion, a box with a square cross-section, 10, is filled with immersion liquid. The video camera, 11, has a frame rate of up to 400 fps. Bubble sizes are calculated using an image-processing procedure in homemade code. The main steps of the algorithm are as follows: 1. The bubble area is chosen on the original image; 2. The image is binarized using an adaptive threshold; 3. Spherical or ellipsoidal particles are sought using a Hough transform. The image processing steps are shown in Figure 3. Figure 3a shows an original image of the bubble flow. The result of image binarization is shown in Figure 3b. The final image with the recognized images of bubbles indicated by blue lines is shown in Figure 3c. Not all of the bubbles are recognized, especially those from the clusters of bubbles, but all groups of bubbles were recognized and we obtain enough statistics for all sizes of bubbles.
For nearly spherical bubbles, the measurement uncertainties depends upon the square of One thousand pairs of PIV images were obtained for both of two cases. One of these cases was a single-phase flow with Re H = 2.09 × 10 4 . In another case, a gas phase with β = 3.5% was added to the flow. Image processing was carried out using the "ActualFlow" software. At first a procedure of subtraction of the mean intensity field averaged over the whole sample range was applied for the raw data. The image processing was carried out by means of iterative cross-correlation algorithm. The overlap of interrogation windows was 75%. The size of the interrogation window was 64 × 64 pixels. Following validation procedures were used: peak validation with the threshold of 1.2, adaptive median 3 × 3 filter and range validation taking into account of maximum velocity on the axis of the small pipe multified by a coefficient 1.2. The uncertainty of the liquid velocity measurements was about 5%.

Bubble Size Measurements
Measurements of the sizes of bubbles are performed using video with shadow illumination prior to the location of the expansion. To reduce the optical distortion, a box with a square cross-section, 10, is filled with immersion liquid. The video camera, 11, has a frame rate of up to 400 fps. Bubble sizes are calculated using an image-processing procedure in homemade code. The main steps of the algorithm are as follows: 1. The bubble area is chosen on the original image; 2. The image is binarized using an adaptive threshold; 3. Spherical or ellipsoidal particles are sought using a Hough transform. The image processing steps are shown in Figure 3. Figure 3a shows an original image of the bubble flow. The result of image binarization is shown in Figure 3b. The final image with the recognized images of bubbles indicated by blue lines is shown in Figure 3c. Not all of the bubbles are recognized, especially those from the clusters of bubbles, but all groups of bubbles were recognized and we obtain enough statistics for all sizes of bubbles.

Physical Model
The modeling of gas bubbles is treated by the Eulerian approach that considers the particulate phase as a continuous medium with properties analogous to those of a fluid (water) [17][18][19]. The Eulerian approach is based on kinetic equations for a one-point probably density function (PDF) of bubbles coordinates, velocity, and temperature in the turbulent Gaussian fluid flow fields [12,[20][21][22]. The system of equations for the carrier phase (liquid), the turbulence model and the equations for the dispersed phase (gas bubbles) are written in a form appropriate to axisymmetric flow, but for reason of brevity, are given below using the vector analysis form [14]:

The Mean Governing Equations for the Carrier Fluid Phase
The mean fluid (water) flow is treated as a steady-state, incompressible and axisymmetric flow. The set of RANS equations for the carrier fluid phase consists of mass conservation, two-momentum, and energy equations. (1) where Pin denotes the mean pressure at the liquid (water)-gas bubble interface; g is gravitational acceleration; τ, ' ' u u and θ u are the viscous stress, Reynolds stress and turbulent heat flux, respectively; M is the volume interphase interaction term (which will be defined in Section 3.5) and gut is the coefficient of involvement of dispersed phase into the large-eddy temperature fluctuation motion of fluid phase [12]: where TLP,t is the time of interaction of bubbles with thermal turbulent eddies TLP,t = 0.6k/ε [12]. The similar simplifications were adopted in many papers. The non-steady-state pressure on the For nearly spherical bubbles, the measurement uncertainties depends upon the square of projection in the measurement area and upon the pixel resolution. The determination error for the bubble diameters for our cases (1-3 mm) is in the range of 3 ÷ 7% depending on the sizes of bubbles.

Physical Model
The modeling of gas bubbles is treated by the Eulerian approach that considers the particulate phase as a continuous medium with properties analogous to those of a fluid (water) [17][18][19]. The Eulerian approach is based on kinetic equations for a one-point probably density function (PDF) of bubbles coordinates, velocity, and temperature in the turbulent Gaussian fluid flow fields [12,[20][21][22]. The system of equations for the carrier phase (liquid), the turbulence model and the equations for the dispersed phase (gas bubbles) are written in a form appropriate to axisymmetric flow, but for reason of brevity, are given below using the vector analysis form [14]:

The Mean Governing Equations for the Carrier Fluid Phase
The mean fluid (water) flow is treated as a steady-state, incompressible and axisymmetric flow. The set of RANS equations for the carrier fluid phase consists of mass conservation, two-momentum, and energy equations.
where P in denotes the mean pressure at the liquid (water)-gas bubble interface; g is gravitational acceleration; τ, u u and uθ are the viscous stress, Reynolds stress and turbulent heat flux, respectively; M l = −M b is the volume interphase interaction term (which will be defined in Section 3.5) Energies 2019, 12, 2735 6 of 18 and g ut is the coefficient of involvement of dispersed phase into the large-eddy temperature fluctuation motion of fluid phase [12]: where T LP,t is the time of interaction of bubbles with thermal turbulent eddies T LP,t = 0.6k/ε [12]. The similar simplifications were adopted in many papers. The non-steady-state pressure on the gas "liquid-bubble" interface is affected by the turbulence of the fluid phase and a non-stationary interfacial slip velocity [12]. The influence of turbulence by now is rather difficult to calculate. A simpler problem in the mathematical plan is the calculation of the effect of the averaged interphase velocity [23][24][25][26][27]. When calculating the pressure on the interphase surface, one can use the Lamb' expression for the potential flow of a particle by a liquid flow [28]: In this study, the value of the constant is assumed C = 0.5 [25,26]. The pressure value at the "liquid-bubble" interface is equal to that one inside the gas bubble for P b << P, P in = P b by [23,24]. Pr T = 0.85 is the turbulent Prandtl number. The Boussinesq hypothesis is used for calculation of turbulent heat flux in the fluid phase

The Turbulence Model of the Carrier Phase (Liquid)
In the present study, the low-Reynolds number elliptic blending second-moment closure from [29] was employed: Here, T T and L T are the time of turbulent macro-scale and turbulent macro-scale length, δ is the Kronecker symbol, χ is the blending coefficient. The value χ changes from zero at the wall to unity far from it.
The terms S k and S ε determine the turbulence generation in the carrier fluid (liquid) by the bubbles presence [23]: where C 3 = 1.92, C 4 = 0.1, C ε4 = 1.44, C D , α and d are the drag coefficient, the void fraction and the gas bubble diameter. The turbulent kinetic energy of the carrier fluid of the two-phase bubbly flow k is calculated using the formula: k = 0.5 u u .

The System of Basic Equations for the Dispersed Phase (Gas Bubbles)
The set of mean governing equations for the dispersed phase has the form [14]: Here, t is the time, are the turbulent diffusion tensor and turbulent heat flux in the dispersed phase [20]; are dynamic and thermal relaxation times [20], the added mass coefficient is C VM = 0.5 [30], and C D is the drag coefficient [31].

Interface Forces
The interfacial force is usually consists of a few components: drag, virtual mass, gravity, lift, turbulent dispersion and wall lubrication: The interphase interaction M l is considered by taking into account the effect of following forces: drag F Drag , added mass F VM , gravity F GA , lift F L , turbulent diffusion (dispersion) F TD and wall lubrication F WL . The interfacial interaction M l in Equation (9) takes into account only the averaged fluid parameters in the paper.
The drag coefficient C D of the bubbles is calculated by the formula [31]: It was obtained the C D in tap water is the same as for solid particles (see paper [14]). The well-known lift force formulation for two-phase flows, which has a positive coefficient C L , acts in the direction of decreasing liquid velocity. The expression for prediction of lift force has the form [32] and C L is the lift coefficient [32]: here f (Eo b ) is the correction function taking into account the bubble deformation, and d H is the the maximum horizontal dimension of the bubble by [33] The coefficient C L changes its sign at a bubble diameter of d = 5.8 mm for our conditions. It was previously shown that this expression of the lift force can be sufficiently used for predictions of bubbly flows in complex geometrical conditions, as example in fuel columns with grid spacers [34]. C TD = 0.1 is the coefficient of turbulent diffusion [23]. In order to handle this behavior of bubbles near the wall an additional wall lubrication force is introduced [35]. The formula for F WL has been modified for the flow in a pipe [36] and the coefficient C W has the form [36]:

Numerical Procedures, Boundary Conditions
The mean transport equations for fluid (liquid) and dispersed (gas bubbles) phases and the turbulence model are solved using a control volumes method on a staggered grid. The QUICK scheme is used to approximate the convective transport. The central difference scheme of the second-order accuracy is performed for the diffusion terms. The SIMPLEC algorithm is employed for coupling velocity and pressure.
All velocity components, temperatures of the phases and turbulence levels are uniform at the inlet. The symmetry conditions are set on the pipe axis for gas and dispersed phases. No-slip conditions are set on the wall surface for the carrier phase. At the outlet edge, the computational domain condition ∂φ/∂r = 0 is set for all variables. The first cell was located at a distance y + = yu*/υ = 0.3-0.5 for all simulations, where u* is the friction velocity obtained for the one-phase flow in the inlet. At least 10 control volumes were generated to ensure resolution of the mean velocity field and turbulence quantities in the viscosity-affected near-wall region (y + < 10).
The correctness of a numerical simulation is strongly dependent upon the quality of the grid. Grid sensitivity study was performed to determine the optimum grid resolution to give the mesh-independent solution. Grid dependence was verified for three different grid sizes: 128 × 50, 256 × 100, 450 × 200, and 550 × 200 control volumes (CVs) in the axial and radial directions (see Figure 4). The distributions of the Nusselt number along the longitudinal coordinate are presented in Figure 4a. The Nusselt number was calculated according to the following relationship for the case q W = const: here H is the step height, q W is the heat flux density for the pipe wall, λ is the fluid's coefficient of heat conductivity, and T m2 and T m1 are the mean liquid temperatures at the pipe axis in the outlet and inlet sections respectively, while the mean liquid temperatures are calculated using the formula: The results obtained with the grid of 128 × 50 cells deviate from those using the other two grids for all pipe lengths.    The results obtained with the grid of 128 × 50 cells deviate from those using the other two grids for all pipe lengths. The predicted values of the Nusselt number obtained for grids with 256 × 100, 450 × 200, and 550 × 200 control volumes (CVs) almost overlap each other at all pipe lengths. The basic grid with 256 × 100 CVs was chosen for all numerical investigations performed. The profiles of turbulent kinetic energy (TKE) in radial direction are shown in Figure 4b

Results and Discussion
We estimated the effect of bubble break up and coalescence on heat transfer modification and we obtained that the influence is not significant (up to 10%) for the conditions of the study. Therefore, we have not been taking into account these effects and it allows us simplifying the mathematical model. The initial size of spherical bubbles is obtained from our measurements.

Flow Structure
The measured and predicted profiles for the mean velocities of the carrier fluid phase and gas bubbles are presented in Figure 6 at a few stations downstream of the pipe's abrupt expansion. The first two sections are situated in the flow recirculation region and two other stations are set in the flow relaxation zone. Open symbols are the one-phase fluid flow and solid symbols are the carrier fluid in the bubbly flow. Solid lines are the predicted data for the fluid phase and dashed curves are the predicted results for the air bubbles. The liquid velocity for bubbly flows (2, 4) is higher than those for the one-phase fluid flow (1, 3). The predicted mean axial velocities for air bubbles (5) are higher than those ones for the liquid velocity in the presence of air bubbles due to the upward direction of the two-phase flow. Thus, it is seen that the axial velocities of the liquid in bubbly flow have negative values at first two sections (x/H = 4 and 8). The carrier fluid in the bubbly flow forms a zone with negative magnitude of mean axial velocity in the near-wall area. These conclusions agree with those for the one-phase separated flow [1,2,37].
We estimated the effect of bubble break up and coalescence on heat transfer modification and we obtained that the influence is not significant (up to 10%) for the conditions of the study. Therefore, we have not been taking into account these effects and it allows us simplifying the mathematical model. The initial size of spherical bubbles is obtained from our measurements.

Flow Structure
The measured and predicted profiles for the mean velocities of the carrier fluid phase and gas bubbles are presented in Figure 6 at a few stations downstream of the pipe's abrupt expansion. The first two sections are situated in the flow recirculation region and two other stations are set in the flow relaxation zone. Open symbols are the one-phase fluid flow and solid symbols are the carrier fluid in the bubbly flow. Solid lines are the predicted data for the fluid phase and dashed curves are the predicted results for the air bubbles. The liquid velocity for bubbly flows (2,4) is higher than those for the one-phase fluid flow (1, 3). The predicted mean axial velocities for air bubbles (5) are higher than those ones for the liquid velocity in the presence of air bubbles due to the upward direction of the two-phase flow. Thus, it is seen that the axial velocities of the liquid in bubbly flow have negative values at first two sections (x/H = 4 and 8). The carrier fluid in the bubbly flow forms a zone with negative magnitude of mean axial velocity in the near-wall area. These conclusions agree with those for the one-phase separated flow [1,2,37].   Figure 7a shows the normalized kinetic energy of turbulence for the fluid in the two-phase flow profiles. The unpredicted tangential normal stress w'w' was equal to the radial normal stress v'v': A significant increase in the fluid carrier phase turbulent kinetic energy (TKE) is observed in the case of bubbly flow. The increase in the intensity of the fluid phase turbulence level in the two-phase flow (up to 20-30%) is shown in comparison with the single-phase flow (1). The additional production of carrier fluid phase turbulence is explained by vortex formation upon streamlining of the gas bubbles by the carrier fluid flow. As is expected, the fluid's maximum turbulent kinetic energy observes in the shear mixing layer. The profiles of TKE at a small value of β agree qualitatively with those for the case of one-phase flow in a pipe with abrupt expansion.  Figure 7a shows the normalized kinetic energy of turbulence for the fluid in the two-phase flow profiles. The unpredicted tangential normal stress w'w' was equal to the radial normal stress v'v': A significant increase in the fluid carrier phase turbulent kinetic energy (TKE) is observed in the case of bubbly flow. The increase in the intensity of the fluid phase turbulence level in the two-phase flow (up to 20-30%) is shown in comparison with the single-phase flow (1). The additional production of carrier fluid phase turbulence is explained by vortex formation upon streamlining of the gas bubbles by the carrier fluid flow. As is expected, the fluid's maximum turbulent kinetic energy observes in the shear mixing layer. The profiles of TKE at a small value of β agree qualitatively with those for the case of one-phase flow in a pipe with abrupt expansion. The distributions of the wall friction coefficient, along the pipe length are shown in Figure 7b. The increase in the β leads to a significant increase in the absolute value of the wall friction coefficient. Note that the minimum value of friction on the wall, located in the recirculation area, and it shifts slightly toward the inlet cross-section with increasing concentration.  The distributions of the wall friction coefficient,

Heat Transfer
along the pipe length are shown in Figure 7b. The increase in the β leads to a significant increase in the absolute value of the wall friction coefficient. Note that the minimum value of friction on the wall, located in the recirculation area, and it shifts slightly toward the inlet cross-section with increasing concentration.  The increase in heat transfer reaches almost three times for the case of ReH = 1.02 × 10 4 and β = 10.1% according to the measurements and simulations. The increased heat transfer coefficient is caused by significant deformation in the velocity profiles of the two-phase flow compared to the single-phase flow and by the increase in the liquid velocity gradient near the pipe wall (see Figure 7). A significant growth in wall friction is obtained in the bubbly flow, even when very small gas concentrations are present, as previously detected for turbulent [37] and laminar [38] flow regimes. The heat transfer enhancement is numerically predicted mainly in the flow relaxation zone. Actually, these conclusions are in agreement with the measured one. Only a few bubbles penetrate into the flow recirculation region and bubbles are available only in the core zone and shear layer region. We obtain a small effect of the gas bubbles on the measured and predicted heat transfer in the flow recirculation region. Authors of recent numerical work [14] showed that small bubbles (d < 1.5 mm) caused heat transfer intensification over the entire length of the recirculation zone, while the larger ones caused intensification mostly in the flow relaxation region. The bubbles migrate towards the wall after the reattachment point of the flow and accumulate there due to the action of the lift force; it leads to an increase in fluid phase turbulence and heat transfer in this region. The length of the zone of the heat transfer intensification is limited for the case of one-phase fluid flow. It is shown, that relatively small amount of the gas bubbles leads to the significant increase of the region where the intensification of heat transfer can be obtained.

The Comparison with Results of Other Papers
The modeling of bubble distributions across the pipe radius is crucial point for the predictions of interfacial term in turbulent bubbly flows. We did not measure the radial bubble distributions and therefore the developed model was additionally validated against measured and numerical results for isothermal, upward bubbly flows downstream from a sudden pipe expansion. Points are the measurements [11], solid lines represent the LES computations [10], and the dashed lines are authors' simulations. The bulk velocity was Um1 = 1.78 m/s, which corresponds to the Reynolds number ReH = 1.11 × 10 5 . The pipe diameter, before the abrupt, was 2R1 = 50 mm, while after expansion, 2R2 = 100 mm, which corresponded to the step height of H = 25 mm, and ER = (R2/R1) 2 = 4 and the inlet bubble diameter was d = 2 mm. The measurements were carried out at five stations located downstream from the expansion edge x = 70 mm (x/H = 2.8), 130 mm (5.2), 250 mm (10) and 320 mm (12.8). The measurements [11] were performed without interphase heat transfer and realized for very different values of Reynolds number and flow geometry than that one of [15]. The increase in heat transfer reaches almost three times for the case of Re H = 1.02 × 10 4 and β = 10.1% according to the measurements and simulations. The increased heat transfer coefficient is caused by significant deformation in the velocity profiles of the two-phase flow compared to the single-phase flow and by the increase in the liquid velocity gradient near the pipe wall (see Figure 7). A significant growth in wall friction is obtained in the bubbly flow, even when very small gas concentrations are present, as previously detected for turbulent [37] and laminar [38] flow regimes. The heat transfer enhancement is numerically predicted mainly in the flow relaxation zone. Actually, these conclusions are in agreement with the measured one. Only a few bubbles penetrate into the flow recirculation region and bubbles are available only in the core zone and shear layer region. We obtain a small effect of the gas bubbles on the measured and predicted heat transfer in the flow recirculation region. Authors of recent numerical work [14] showed that small bubbles (d < 1.5 mm) caused heat transfer intensification over the entire length of the recirculation zone, while the larger ones caused intensification mostly in the flow relaxation region. The bubbles migrate towards the wall after the reattachment point of the flow and accumulate there due to the action of the lift force; it leads to an increase in fluid phase turbulence and heat transfer in this region. The length of the zone of the heat transfer intensification is limited for the case of one-phase fluid flow. It is shown, that relatively small amount of the gas bubbles leads to the significant increase of the region where the intensification of heat transfer can be obtained.

The Comparison with Results of Other Papers
The modeling of bubble distributions across the pipe radius is crucial point for the predictions of interfacial term in turbulent bubbly flows. We did not measure the radial bubble distributions and therefore the developed model was additionally validated against measured and numerical results for isothermal, upward bubbly flows downstream from a sudden pipe expansion. Points are the measurements [11], solid lines represent the LES computations [10], and the dashed lines are authors' simulations. The bulk velocity was U m1 = 1.78 m/s, which corresponds to the Reynolds number Re H = 1.11 × 10 5 . The pipe diameter, before the abrupt, was 2R 1 = 50 mm, while after expansion, 2R 2 = 100 mm, which corresponded to the step height of H = 25 mm, and ER = (R 2 /R 1 ) 2 = 4 and the inlet bubble diameter was d = 2 mm. The measurements were carried out at five stations located downstream from the expansion edge x = 70 mm (x/H = 2.8), 130 mm (5.2), 250 mm (10) and 320 mm (12.8). The measurements [11] were performed without interphase heat transfer and realized for very different values of Reynolds number and flow geometry than that one of [15]. The development of the axial fluid phase velocity along different sections of the pipe is shown in Figure 9a. The maximal difference between measurements of [11], LES results of [10] and our RANS predictions are observed in the flow recirculation region (the maximal discrepancy is up to 15%). The predicted radial profiles for the liquid phase axial mean fluid phase velocity (a), liquid velocity fluctuations (b) and bubble volume fraction (c) are presented in Figure 9. The gas bubble volume fraction showed a strongly pronounced maximum in the shear layer for first two sections. It can be explained by turbophoresis (turbulent transport), turbulent diffusion forces [10,22] and bubbles accumulating in zones with a high value for fluid flow pulsations. Figure 9b shows the distributions of axial fluid (liquid) phase velocity pulsations at a few stations along the pipe length. After the bubbly flow reattachment at x/H ≈ 9, along with its relaxation, there is a transition to the turbulent fully developed flow. The magnitude of axial fluid velocity fluctuations in this zone has the highest value too and the shear layer is a "trap" for them. The profiles for the void fractions in the recirculation zone are characterized by an almost zero value of the bubble concentration at first two stations (x/H < 10). It causes heat transfer enhancement predicted in our simulations, lesser than that one in the flow relaxation area. The radial profile for the bubble concentration became more uniform in the region after the flow reattachment (x/H > 10) compared to one in the first two stations (in the region of flow recirculation). The crest of the fluid flow pulsations disappeared and lead to the redistribution of bubbles across the pipe's cross section as they moved toward the wall (see Figure 9c). The values of bubble volume fractions have almost zero magnitude in pipe's near-wall region. The maximal difference between measured and numerical results is up to 20%. The authors results agree well with experiments and LES data. The development of the axial fluid phase velocity along different sections of the pipe is shown in Figure 9a. The maximal difference between measurements of [11], LES results of [10] and our RANS predictions are observed in the flow recirculation region (the maximal discrepancy is up to 15%). The predicted radial profiles for the liquid phase axial mean fluid phase velocity (a), liquid velocity fluctuations (b) and bubble volume fraction (c) are presented in Figure 9. The gas bubble volume fraction showed a strongly pronounced maximum in the shear layer for first two sections. It can be explained by turbophoresis (turbulent transport), turbulent diffusion forces [10,22] and bubbles accumulating in zones with a high value for fluid flow pulsations. Figure 9b shows the distributions of axial fluid (liquid) phase velocity pulsations at a few stations along the pipe length. After the bubbly flow reattachment at x/H ≈ 9, along with its relaxation, there is a transition to the turbulent fully developed flow. The magnitude of axial fluid velocity fluctuations in this zone has the highest value too and the shear layer is a "trap" for them. The profiles for the void fractions in the recirculation zone are characterized by an almost zero value of the bubble concentration at first two stations (x/H < 10). It causes heat transfer enhancement predicted in our simulations, lesser than that one in the flow relaxation area. The radial profile for the bubble concentration became more uniform in the region after the flow reattachment (x/H > 10) compared to one in the first two stations (in the region of flow recirculation). The crest of the fluid flow pulsations disappeared and lead to the redistribution of bubbles across the pipe's cross section as they moved toward the wall (see Figure  9c). The values of bubble volume fractions have almost zero magnitude in pipe's near-wall region. The maximal difference between measured and numerical results is up to 20%. The authors results agree well with experiments and LES data.

Conclusions
The turbulent flow patterns and the heat transfer for a bubbly turbulent flow in a pipe with sudden expansion were experimentally and numerically investigated. An experimental study of structure of bubbly flow was carried out using PIV-LIF approach and shadow photography. Experimental and numerical investigations were performed using a range of Reynolds numbers of ReH = (1-2) × 10 4 and β = 0−10%. The numerical model was employed the Eulerian approach. The set of RANS equations was used for modelling two-phase bubbly flows. The turbulence of the carrier liquid phase was predicted using the Reynolds stress model. A significant growth (up to 30%) in wall friction in bubbly flow downstream of pipe with abrupt expansion is obtained. The mean and fluctuating flow structures in bubbly flow with small values of β ≤ 10% is similar to that of a one-phase fluid flow. The increase in the intensity of the fluid phase turbulence level in the two-phase flow (up to 20-30%) is shown. The peak of axial and radial fluctuations of the carrier fluid (liquid) velocity in the bubbly flow is observed in the shear layer. It is experimentally and numerically shown the heat transfer in bubbly flow is significantly enhances (up to three times). This effect is augmented with increasing gas volumetric flow rate ratios β. The main enhancement in heat transfer is observed after the point of flow reattachment. In general, the Nusselt number distributions in the bubbly flow have qualitatively similar character as that one for singlephase fluid turbulent separated flows.
Two regions can be found for the influence of bubbles on the flow. Similar behaviour of twophase flow was previously found by the authors in the model of pressurized water reactor downstream a spacer grid [39]. In the recirculation zone the heat transfer is determined by large vortices structures. The effect of bubbles on the heat transfer in this region is minimal. The main increase in heat transfer is observed in the flow relaxation region after the point of flow reattachment. It is experimentally and numerically shown that the addition of air bubbles causes a significant increase in the heat transfer rate (up to three times), and these effects increase with increasing gas volumetric flow rate ratios. Intensification of the heat transfer in this region relates to the additional flow turbulisation by bubbles and the reorganisation of the liquid velocity profile which leads to the increase of the velocity gradient in the near wall region. Early a similar degree of the increasing of wall shear stress and heat transfer coefficient in turbulent two-phase bubbly flows was found in [40,41].

Conclusions
The turbulent flow patterns and the heat transfer for a bubbly turbulent flow in a pipe with sudden expansion were experimentally and numerically investigated. An experimental study of structure of bubbly flow was carried out using PIV-LIF approach and shadow photography. Experimental and numerical investigations were performed using a range of Reynolds numbers of Re H = (1-2) × 10 4 and β = 0−10%. The numerical model was employed the Eulerian approach. The set of RANS equations was used for modelling two-phase bubbly flows. The turbulence of the carrier liquid phase was predicted using the Reynolds stress model. A significant growth (up to 30%) in wall friction in bubbly flow downstream of pipe with abrupt expansion is obtained. The mean and fluctuating flow structures in bubbly flow with small values of β ≤ 10% is similar to that of a one-phase fluid flow. The increase in the intensity of the fluid phase turbulence level in the two-phase flow (up to 20-30%) is shown. The peak of axial and radial fluctuations of the carrier fluid (liquid) velocity in the bubbly flow is observed in the shear layer. It is experimentally and numerically shown the heat transfer in bubbly flow is significantly enhances (up to three times). This effect is augmented with increasing gas volumetric flow rate ratios β. The main enhancement in heat transfer is observed after the point of flow reattachment. In general, the Nusselt number distributions in the bubbly flow have qualitatively similar character as that one for single-phase fluid turbulent separated flows.
Two regions can be found for the influence of bubbles on the flow. Similar behaviour of two-phase flow was previously found by the authors in the model of pressurized water reactor downstream a spacer grid [39]. In the recirculation zone the heat transfer is determined by large vortices structures. The effect of bubbles on the heat transfer in this region is minimal. The main increase in heat transfer is observed in the flow relaxation region after the point of flow reattachment. It is experimentally and numerically shown that the addition of air bubbles causes a significant increase in the heat transfer rate (up to three times), and these effects increase with increasing gas volumetric flow rate ratios. Intensification of the heat transfer in this region relates to the additional flow turbulisation by bubbles and the reorganisation of the liquid velocity profile which leads to the increase of the velocity gradient in the near wall region. Early a similar degree of the increasing of wall shear stress and heat transfer coefficient in turbulent two-phase bubbly flows was found in [40,41].