**1. 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–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 regime. Gas concentration increased immediately after the cross-section of the flow separation. Measurements by [7] were performed in an upward, vertical flow of a liquid (water) and CO2 bubbles. The distribution of the mean carrier fluid velocity was strongly dependent on the bubbles' diameter. A review of the literature for two-phase flow in sudden expansions and contractions was performed in [8]. Wang et al. [8] have found that none of the existing eight correlations can accurately predict the experimental database. Most of these correlations greatly overpredicted the results for mini test sections. Some of the correlations also underpredicted the data for large test sections.

There are only a few papers that consider the numerical modeling of isothermal bubbly flows in a pipe with sudden expansion [9,10]. Krepper et al. [9] presented a mathematical model for polydisperse flow in a vertical pipe with sudden expansion. The numerical predictions for bubbly pipe flow with an obstacle demonstrated very difficult nature of such flows [9]. The Eulerian model was developed to simulate bubbly flows in a pipe with sudden expansion in [10]. Predictions were carried out using the CFD commercial package STAR-CCM+. The developed model was validated against measurements of Bel Fdhila [11]. The eddy-viscosity k−<sup>ε</sup> model showed poor reproduction of the bubbly flow structure because it disregarded much of the information on fluid flow turbulence. The Second-Moment Closure and LES methods captured the main characteristic turbulent length-scales responsible for the diffusion of the bubble gas phase, and moreover, these methods predicted the flow patterns and bubble distributions across the pipe section [10].

Recently published numerical works [12,13] investigated heat transfer in bubbly flows without an abrupt pipe or channel expansion. These studies showed the significant influence of bubble diameter and gas volumetric flow rate ratios on the flow patterns and heat transfer. Heat transfer in bubbly flow increased up to three times that of one-phase fluid flow. The liquid turbulence was modeled using an isotropic k−<sup>ε</sup> model [12]. The DNS of bubbly flow with heat transfer was performed [13]. The numerical [14,15] and experimental papers [15] concerned study of bubbly flows with a sudden pipe expansion. The numerical model [14] employed the Eulerian approach. The turbulence of the carrier fluid flow was predicted with the use of a Reynolds stress model. The effects of inlet gas concentrations, Reynolds numbers, and bubble diameters on the flow structure and heat transfer were investigated. The measurements of heat transfer in gas–liquid flow in the pipe with sudden expansion and validation analysis of model [14] was performed [15].

The aim of this work is to carry out a further validation of the numerical model [14] for predictions of the flow structure in pipe sudden expansion. Authors of [15] used for comparison analysis only their measured results of heat transfer in the gas–liquid flow downstream of the pipe with sudden expansion. In the paper we performed the validations with our experimental results of mean flow structure, measurements of [11] and LES simulations of [10]. The measured and predicted results for the mean axial carrier fluid (liquid) and gas (air bubbles) phases in a few stations along the pipe length are presented. The LES predictions of [10] and measurements of [11] were performed for much larger Reynolds number (Re = 1.1 × 105) and another pipe geometry (2*R*1 = 50 mm, 2*R*2 = 100 mm, *H* = 25 mm, and ER = (*R*2/*R*1)<sup>2</sup> = 4). It has expanded the possibility of using of the model. We have also presented the numerical results of the flow structure predictions (turbulent kinetic energy profiles and wall friction coefficient distributions for various β). All these abovementioned numerical results did not present in [15].

This study may be of interest to scientists and engineers dealing with the enhancement of heat transfer in power equipment. The paper does not have any direct practical application. The article shows the physical mechanisms of the control on flow structure and heat transfer enhancement in bubbly flow in a pipe sudden expansion.

#### **2. Experimental Facility and Image Processing Methods**

#### *2.1. Experimental Facility*

The experimental facility consists of a closed circuit in a liquid and an open circuit in the gas phase (see Figure 1). The test liquid, that is, distilled water, is pumped from the main tank, 1, to the test Section 2, through the water supply line, 3. The volume of liquid in the main tank is 40 L. To 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 managemen<sup>t</sup> 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—thepipeoflargediameter(afterexpansion);17—downcomer.

A Grundfos CHI 4–40 (Ql ≤ 4.5 m<sup>3</sup>/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)<sup>2</sup> = 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.

#### *2.2. 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 di fferent 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.

**Figure 2.** The photographs of bubbles in the sudden pipe expansion two-phase flow. Re *H* = 2.09 × 104, 2*R*1 = 15 mm, *H* = 13.5 mm, 2 *R*2 = 42 mm, ER = (*R*2/*R*1) 2 = 7.84. (**a**)—β = 0.6%, (**b**)—1.8, (**c**)—3.5; (**d**)—5.2.

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 × 104. 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 coe fficient 1.2. The uncertainty of the liquid velocity measurements was about 5%.

#### *2.3. 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.

**Figure 3.** Example of picture processing: (**a**)—original picture; (**b**)—binared picture; (**c**)—picture after determination of bubbles.

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.

#### **3. Mathematical Model for Two-Phase Bubbly Flow**

#### *3.1. 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–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–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]:

#### *3.2. 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.

> τ −

$$\begin{split} \nabla \cdot (a\_l \rho \mathbf{U}) &= 0, \\ \nabla \cdot (a\_l \rho \mathbf{U} \mathbf{U}) &= a\_l (-\nabla P + \rho \mathbf{g}) + \nabla \cdot (a\_l \mathbf{r}) - \nabla \cdot (a\_l \rho (\mathbf{u'} \mathbf{u'})) + (P - P\_{\rm in}) \nabla \cdot a\_l + \mathbf{M}\_l, \\ \nabla \cdot (a\_l \rho \mathbf{C}\_P \mathbf{U} \mathbf{T}) &= \nabla \cdot [a\_l \nabla (\lambda \mathbf{T})] - \nabla \cdot (a\_l \rho \mathbf{C}\_P \langle \mathbf{u} \theta \rangle) + \frac{6 \mathbf{b} (T - T\_b)}{d} + \\ &+ \mathbf{C}\_{\rm Pb} \rho\_b \mathbf{g}\_{\rm at} (\mathbf{u} \theta) \nabla \cdot a\_l \end{split} \tag{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; *Ml* = −*M<sup>b</sup>* 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]:

$$g\_{\rm ut} = T\_{\rm LP,t}/\tau\_{\rm \Theta} - 1 + \exp(T\_{\rm LP,t}/\tau\_{\rm \Theta}),\\\tau\_{\rm \Theta} = C\_{\rm Pb} \rho\_b d^2 / (12 \,\text{Y} \lambda),\\Y = 1 + 0.3 \text{Re}\_b^{1/2} \text{Pr}^{1/3} \tag{2}$$

where *TLP,t* is the time of interaction of bubbles with thermal turbulent eddies *TLP,t* = 0.6*k*/ε [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–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]:

$$P\_{\rm int} = P\_b = P - C\rho\alpha |\mathbf{U}\_R|^2. \tag{3}$$

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 *Pb* << *P*, *Pin* = *Pb* by [23,24]. Pr<sup>T</sup> = 0.85 is the turbulent Prandtl number. The Boussinesq hypothesis is used for calculation of turbulent heat flux in the fluid phase

$$
\left\langle \mathbf{u}\_{\mathrm{j}} \boldsymbol{\Theta} \right\rangle = -\frac{\mathbf{v}^{\mathrm{T}}}{\mathrm{Pr}^{\mathrm{T}}} \frac{\partial \mathbf{T}}{\partial \mathbf{x}\_{\mathrm{j}}}.
$$

#### *3.3. 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:

$$\begin{split} \nabla(\mathbf{u}|\rho(\mathbf{u}'\mathbf{u}')) &= \mathbf{a}[\mathbf{P} + \Phi - \varepsilon\_{\rm d} + \nabla\Big(\mathbf{v}\delta + \frac{\mathbf{C}\_{\rm u}\mathbf{T}^{\rm T}}{\sigma\_{\rm k}}\langle\mathbf{u}'\mathbf{u}'\rangle\Big)\nabla\langle\mathbf{u}'\mathbf{u}'\rangle] + \mathbf{S}\_{\rm k} \\ \nabla(\mathbf{a}\_{\rm l}p\varepsilon) &= \mathbf{a}[\Big[\frac{1}{\mathbf{T}^{\rm T}}(\mathbf{C}\_{\varepsilon1}\mathbf{P}\_{2} - \mathbf{C}\_{\varepsilon2}\varepsilon) + \nabla\Big(\mathbf{v}\delta + \frac{\mathbf{C}\_{\rm u}\mathbf{T}^{\rm T}}{\sigma\_{\rm k}}\langle\mathbf{u}'\mathbf{u}'\rangle\Big)\nabla\varepsilon\Big] + \mathbf{S}\_{\varepsilon} \\ \chi - (\mathbf{L}\_{\Gamma})^{2}\nabla^{2}\chi &= 1. \end{split} \tag{4}$$

Here, T<sup>T</sup> and LT 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.

$$\mathbf{T}^{\mathrm{T}} = \max\{\mathbf{k}/\varepsilon; 6\sqrt{\mathbf{v}/\varepsilon}\},\\L\_{\mathrm{T}} = 0.45\left[\frac{k^{3/2}}{\varepsilon}; 80\left(\frac{\nu}{\varepsilon}\right)^{3/4}\right] \tag{5}$$

The terms Sk and Sε determine the turbulence generation in the carrier fluid (liquid) by the bubbles presence [23]:

$$\mathbf{S}\_{\mathbf{k}} = \mathbf{C}\_{4} \frac{3}{4} \mathbf{C}\_{\mathbf{D}} \alpha \frac{\rho \left| \mathbf{U}\_{\mathbf{R}} \right|^{3}}{\mathbf{d}}, \mathbf{S}\_{\varepsilon} = \mathbf{C}\_{3} \frac{\mathbf{k}}{\varepsilon} \mathbf{S}\_{\mathbf{k}\prime} \tag{6}$$

where C3 = 1.92, C4 = 0.1, Cε*<sup>4</sup>* = 1.44, CD, α 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 .

#### *3.4. 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]:

$$\begin{split} \frac{\frac{\partial (\mathbf{a}\_{\rm b} \mathbf{b}\_{\rm b} \mathbf{U}\_{\rm b})}{\partial \mathbf{t}}}{\frac{\partial (\mathbf{a}\_{\rm b} \mathbf{b}\_{\rm b} \mathbf{U}\_{\rm b})}{\partial \mathbf{t}}} &= -\nabla (\mathbf{a}\_{\rm b} \mathbf{P}\_{\rm b}) - \nabla (\mathbf{p}\_{\rm b} \mathbf{E} \langle \mathbf{u}' \mathbf{u}' \rangle) + \mathbf{M}^{\rm b}, \\\frac{\frac{\mathbf{D} (\mathbf{a}\_{\rm b} \mathbf{p}\_{\rm b} \mathbf{C}\_{\rm b} \mathbf{U}\_{\rm b} \mathbf{T}\_{\rm b})}{\mathbf{D}t}}{\mathbf{D}t} &= -\mathbf{h} \mathbf{c}\_{\rm b} (\mathbf{T} - \mathbf{T}\_{\rm b}) \frac{\rho\_{\rm b}}{\tau\_{\rm 0}} + (1 - \alpha\_{\rm b}) \rho \mathbf{C}\_{\rm P} \frac{\mathbf{D} \mathbf{T}}{\mathbf{D}t} - \frac{\mathbf{D}\_{\rm b}^{\rm B}}{\tau\_{\rm 0}} \nabla \alpha\_{\rm b}, \\\rho\_{\rm b} &= \mathbf{P}\_{\rm b} / \{\overline{\mathbf{R}}\_{\rm b} \mathbf{T}\_{\rm b}\}. \end{split} \tag{7}$$

Here, *t* is the time, DφDt = ∂φ∂t + **U**∇φ is the material derivative, E is the expression of [20]

$$\mathbf{E} = \frac{(1 - \mathbf{A})(1 - \mathbf{A}\Omega)}{1 + \Omega}, \mathbf{A} = 1.5\mathfrak{p}\_0/(1 + 0.5\mathfrak{p}\_0), \ \mathfrak{p}\_0 = \mathfrak{p}/\mathfrak{p}\_{\mathfrak{p}}$$

Ω = τ/TLP is the dimensionless timescale, Db and *D*Θ*b* are the turbulent diffusion tensor and turbulent heat flux in the dispersed phase [20];

$$
\pi = \frac{4\mathbf{d}(1 + \mathbf{C}\_{\text{VM}}\rho\_0)}{3\mathbf{C}\_{\text{D}}|\mathbf{U}\_{\text{R}}|\rho\_0}, \pi\_{\Theta} = \mathbf{C}\_{\text{Pb}}\rho\_{\text{b}}\mathbf{d}^2/(6\lambda\mathbf{Y})\tag{8}
$$

are dynamic and thermal relaxation times [20], the added mass coefficient is CVM = 0.5 [30], and CD is the drag coefficient [31].

#### *3.5. Interface Forces*

The interfacial force is usually consists of a few components: drag, virtual mass, gravity, lift, turbulent dispersion and wall lubrication:

$$\begin{split} \mathbf{M}\_{l} &= -\mathbf{M}\_{b} = \mathbf{F}\_{D\mathbf{n}\mathbf{g}} + \mathbf{F}\_{VM} + \mathbf{F}\_{\mathbf{GA}} + \mathbf{F}\_{L} + \mathbf{F}\_{TD} + \mathbf{F}\_{\mathbf{WL}}, \\ \mathbf{F}\_{\mathbf{Drag}} &= \frac{\alpha\_{b}\mathbf{U}\_{R}}{\tau}, \mathbf{F}\_{VM} = \alpha\_{b}\mathbf{A}\mathbf{L}\nabla \cdot \mathbf{U}, \mathbf{F}\_{\mathbf{GA}} = \frac{(1-\rho\_{b})\alpha\_{b}}{1+\mathbb{C}\_{VM\theta}\mathbf{0}}\mathbf{g}, \\ \mathbf{F}\_{L} &= \frac{\mathbf{C}\_{L}\alpha\_{b}\rho\_{0}}{1+\mathbb{C}\_{VM\theta}\mathbf{0}}\mathbf{U}\_{R} \times (\nabla \times \mathbf{U}), \mathbf{F}\_{TD} = -\frac{\mathbf{C}\_{TD}\alpha\_{b}}{\mathbf{S}\mathbf{c}^{T}}\frac{\mu^{T}}{\tau} \left(\frac{\nabla \cdot \mathbf{a}\_{b}}{a\_{b}} - \frac{\nabla \cdot \mathbf{a}\_{l}}{a\_{l}}\right), \\ \mathbf{F}\_{\mathbf{WL}} &= \mathbf{C}\_{W}\frac{d\rho\_{0}a\_{b}|\mathbf{U}\_{R}|^{2}}{2} \Big(\frac{1}{y^{2}} - \frac{1}{(2\mathbb{R}-y)^{2}}\Big). \end{split} \tag{9}$$

The interphase interaction *Ml* is considered by taking into account the effect of following forces: drag *<sup>F</sup>Drag*, added mass *FVM*, gravity *FGA*, lift *FL*, turbulent diffusion (dispersion) *FTD* and wall lubrication *FWL*. The interfacial interaction *Ml* in Equation (9) takes into account only the averaged fluid parameters in the paper.

The drag coefficient *CD* of the bubbles is calculated by the formula [31]:

$$\mathbf{C}\_{D} = \begin{pmatrix} \frac{24}{\text{Re}\mathfrak{g}} \left( 1 + \frac{3}{16} \text{Re}\_{b}^{0.687} \right) , \text{ Re}\_{b} \le 10^{3} \\\ 0.44 , \text{ Re}\_{b} > 10^{3} \end{pmatrix} . \tag{10}$$

It was obtained the *CD* 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 *CL*, acts in the direction of decreasing liquid velocity. The expression for prediction of lift force has the form [32] and *CL* is the lift coefficient [32]:

$$\mathbf{C}\_{L} = \begin{cases} \min\left[0.288 \tanh(0.121 \text{Re}\_{b}), f(\text{Eo}\_{b})\right], \text{Eo}\_{b} \le 4\\ \qquad f(\text{Eo}\_{b}), \ 4 \le \text{Eo}\_{b} \le 10\\ \qquad -0.27, \ \text{Eo}\_{b} > 10 \end{cases} \tag{11}$$

here *f*(Eo*b*) is the correction function taking into account the bubble deformation, and *dH* is the the maximum horizontal dimension of the bubble by [33]

$$f(\text{Fo}\_b) = 0.0011\text{Eo}\_b^3 - 0.0159\text{Eo}\_b^2 - 0.0204\text{Eo}\_b + 0.474,\\d\_H = d\left(1 + 0.163\text{Eo}^{0.757}\right)^{1/3} \tag{12}$$

The coefficient *CL* 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]. *CTD* = 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 *CW* has the form [36]:

$$\mathcal{C}\_{W} = \begin{cases} \exp(-0.933\text{Eo} + 0.179), & 1 \le \text{Eo} \le 5\\ 0.007\text{Eo} + 0.04, & 5 \le \text{Eo} \le 33 \end{cases} \tag{13}$$

#### **4. 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 *qW* = const:

$$\mathbf{Nu} = Hq\_{\rm IN} / \left[ \lambda \left( T\_{m2} - T\_{m1} \right) \right] \tag{14}$$

here *H* is the step height, *qW* is the heat flux density for the pipe wall, λ is the fluid's coefficient of heat conductivity, and *Tm*2 and *Tm*1 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:

$$T\_m = \frac{2}{\mathcal{U}\_1 \mathcal{R}\_1^2} \int\_0^{\mathcal{R}\_2} T \mathcal{U} r dr \tag{15}$$

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 for a few grids. The predicted

values of TKE obtained for grids with 256 × 100, 450 × 200, and 550 × 200 control volumes almost overlap each other over the pipe radius.

**Figure 4.** Grid independence tests for local heat transfer along the longitudinal coordinate (**a**) and radial profile of kinetic energy of turbulence at x/H = 5 (**b**). Points (1) are authors' measured data, curves (2–5) are authors' simulations results. ReH = 2.09 × 104, 2R1 = 15 mm, H = 13.5 mm, 2R2 = 42 mm, ER = (R2/R1)<sup>2</sup> = 7.84, x/H = 5, qW = const = 8.95 kW/m2, T1 = Tb1 = 298 K, d = 1.7 mm. 2—128 × 50, 3—256 × 100, 4—450 × 200, 5—550 × 200 CVs.

The solution methodology of the developed model is shown in the Figure 5.

**Figure 5.** The solution methodology of the developed model. 1—Initial values input for all flow parameters, 2—Predictions of all flow parameters (two-fluid model + SMC), 3—Calculation of T and Tb, 4—Prediction of void fraction, 5—Convergence, 6 – Stop.

#### **5. 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.

#### *5.1. 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 6.** The distributions of mean axial fluid and gas bubbles velocities at *x*/*H* = 16 (**a**), 12 (**b**), 8 (**c**) and 4 (**d**). Re *H* = 2.09 × 104, *d* = 1.7 mm, β = 3.5%. Points and curves are authors' measurements and simulations respectively. *1*,*3—*β = 0 (single-phase fluid flow), *2*,*4*—carrier fluid phase (β = 3.5%), *5*—gas bubbles (β = 3.5%).

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':

$$2k = \mathfrak{u}'\mathfrak{u}' + \mathfrak{v}'\mathfrak{v}' + \mathfrak{w}'\mathfrak{w}' \approx \mathfrak{u}'\mathfrak{u}' + 2\mathfrak{v}'\mathfrak{v}'.\tag{16}$$

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 7.** The profiles of TKE of the carrier fluid in radial direction at *x*/*H* = 5 (**a**) and wall friction coefficient (**b**). Re*H* = 2.09 × 104, 2*R*1 = 15 mm, *H* = 13.5 mm, 2*R*2 = 42 mm, *d* = 1.7 mm. *1—*β = 0 (one-phase water flow), *2*—3.5%, *3*—5.2%.

The distributions of the wall friction coefficient,

$$\mathbb{C}\_f = 2\tau\_W / \left(\rho l I\_{m1}^2\right) . \tag{17}$$

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.

#### *5.2. Heat Transfer*

Figure 8 shows distributions of the parameter for heat transfer enhancement Nu/Nu0,max in bubbly flow in the pipe with sudden expansion with various β. The local Nusselt number is calculated using the Formula (14). Here Nu0,max is the maximal value of the Nusselt number for the one-phase flow. Points and lines are the results of measurements and numerical simulations performed by the authors, and the dashed lines represent the calculation for a one-phase flow (β = 0) under otherwise identical conditions. For visual clarity, only every seventh data point in the experimental data is shown in the figure.

**Figure 8.** The distributions of heat transfer enhancement ratios along the axial coordinate at Re*H* = 1.02 × 104. Points and lines are authors' measurements and computations respectively. *d* = 1.7−2 mm, 2*R*1 = 15 mm, *H* = 13.5 mm, 2*R*2 = 42 mm, ER = (*R*2/*R*1)<sup>2</sup> = 7.84, *qW* = const = 8.95 kW/m2, *T*1 = *Tb*1 = 298 K. *1—*β = 0 (one-phase flow), *2*—1.2%, *3*—3.6, *4*—7, *5*—10.1.

The increase in heat transfer reaches almost three times for the case of Re*H* = 1.02 × 10<sup>4</sup> 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 agreemen<sup>t</sup> 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.

#### *5.3. 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 *Um*1 = 1.78 m/s, which corresponds to the Reynolds number Re*H* = 1.11 × 105. The pipe diameter, before the abrupt, was 2*R*1 = 50 mm, while after expansion, 2*R*2 = 100 mm, which corresponded to the step height of *H* = 25 mm, and ER = (*R*2/*R*1)<sup>2</sup> = 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 di fferent sections of the pipe is shown in Figure 9a. The maximal di fference 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 di ffusion 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 di fference between measured and numerical results is up to 20%. The authors results agree well with experiments and LES data.

(**b**) 

**Figure 9.** *Cont*.

**Figure 9.** The profiles of mean axial fluid velocity (**a**), fluid phase axial velocity fluctuations (**b**) and gas volume fraction (**c**) across the radial coordinate in bubbly isothermal upward flow in the pipe with sudden expansion. *1*—measured results of [15], *2*—LES data [10], *2*—authors' predictions. Re*H* = 1.11 × 105, *Um*1 = 1.78 m/s, *d*1 = 2 mm, 2*R*1 = 50 mm, *H* = 25 mm, 2*R*2 = 100 mm, ER = (*R*2/*R*1)<sup>2</sup> = 4, *T*1 = *Tb*1 = 300 K.

## **6. 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<sup>4</sup> 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].

**Author Contributions:** The developing of the experimental circuit, the test section for optical measurements, and the purposed method of two-phase flow study by means of PIV/PLIF, experimental study, experimental data processing and description of experimental facility were performed by P.L. M.P. and V.T. developed the numerical model and performed the numerical part of the paper.

**Funding:** The experimental study of hydrodynamic characteristics of the flow was performed by financial support of the Russian Foundation for Basic Research (Project No. 18-58-45006), heat transfer measurements were performed under the state contract with IT SB RAS (AAAA-A18-118051690120-2) and all simulations of the paper were carried out under the state contract with IT SB RAS (AAAA-A17-117030310010-9).

**Acknowledgments:** P. Lobanov is grateful to O.N. Kashinsky for providing the heat transfer measurement unit.

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