*Article* **The Effect of a Backward-Facing Step on Flow and Heat Transfer in a Polydispersed Upward Bubbly Duct Flow**

**Tatiana V. Bogatko 1, Aleksandr V. Chinak 2, Ilia A. Evdokimenko 2, Dmitriy V. Kulikov 2, Pavel D. Lobanov <sup>2</sup> and Maksim A. Pakhomov 1,\***


**Abstract:** The experimental and numerical results on the flow structure and heat transfer in a bubbly polydispersed upward duct flow in a backward-facing step are presented. Measurements of the carrier fluid phase velocity and gas bubbles motion are carried out using the PIV/PLIF system. The set of RANS equations is used for modeling the two-phase bubbly flow. Turbulence of the carrier fluid phase is predicted using the Reynolds stress model. The effect of bubble addition on the mean and turbulent flow structure is taken into account. The motion and heat transfer in a dispersed phase is modeled using the Eulerian approach taking into account bubble break-up and coalescence. The method of delta-functions is employed for simulation of distributions of polydispersed gas bubbles. Small bubbles are presented over the entire duct cross-section and the larger bubbles mainly observed in the shear mixing layer and flow core. The recirculation length in the two-phase bubbly flow is up to two times shorter than in the single-phase flow. The position of the heat transfer maximum is located after the reattachment point. The effect of the gas volumetric flow rate ratios on the flow patterns and maximal value of heat transfer in the two-phase flow is studied numerically. The addition of air bubbles results in a significant increase in heat transfer (up to 75%).

**Keywords:** turbulent bubbly flow; backward-facing step; PIV/PLIF measurements; RANS modeling; flow structure; heat transfer

#### **1. Introduction**

Upward gas-liquid bubbly flows in ducts are frequently employed in chemical, nuclear, power engineering and other practical applications. These flows are characterized by a strong interfacial momentum, heat and mass transfer between the carrier fluid phase (water) and dispersed phase (gas bubbles). Bubbly flows in a backward-facing step (BFS) are complicated detachment, recirculation and reattachment of a flow, polydispersity, breakup, coalescence and heat transfer [1,2]. The first papers concerned with the study of such flows were experimental studies of two-phase bubbly flows in ducts and pipes behind a backward-facing step [3,4]. Authors measured the pressure drops, bubble diameter and local void fraction and mean and fluctuational phase velocities.

The complexity of the numerical modeling of such flows is associated with the need to take into account the number of multiscale factors: turbulence of the carrier phase, gas bubble-carrier fluid interaction, bubble break-up and coalescence processes. Only a few papers devoted to the experimental investigations of turbulent bubbly flows in a duct or pipe with sudden expansion without heat transfer were found in the literature [5–7].

A theoretical and experimental investigation of the upward bubbly flow with a sudden pipe expansion was carried out in [5]. The distributions along the pipe length and across

**Citation:** Bogatko, T.V.; Chinak, A.V.; Evdokimenko, I.A.; Kulikov, D.V.; Lobanov, P.D.; Pakhomov, M.A. The Effect of a Backward-Facing Step on Flow and Heat Transfer in a Polydispersed Upward Bubbly Duct Flow. *Water* **2021**, *13*, 2318. https:// doi.org/10.3390/w13172318

Academic Editor: Majid Mohammadian

Received: 29 June 2021 Accepted: 20 August 2021 Published: 24 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the pipe radius of mean and fluctuating axial velocities of the gas bubbles, local void fraction, distribution of bubble diameter, pressure drop and wall friction coefficient were studied in this work. The experimental study of a mixture of air and oil in a horizontal BFS flow was performed in a bubbly flow mode, in a transitional regime from bubbly to annular flows and in an annular regime [6]. It is shown that an increase in the gas volumetric flow rate ratio after the two-phase flow detachment cross-section occurs due to the separation of air from the two-phase flow in the recirculation region.

The flow structure of oil-water and kerosene-water flows in a horizontal conduit behind a sudden contraction and expansion has been experimentally studied [7]. Thick core flow can be revealed for high viscous oil in a sudden expansion. It increases the probability of wall fouling by the oil. The experimental values of contraction and expansion coefficients are found to be lower for oil-water as compared to only single-phase fluid water flow for other conditions being identical. Experimental studies of the pressure drop for a two-phase slug flow in a horizontal channel with a sudden expansion were carried out in work [8]. An electrochemical method was used to measure the wall shear stress. It has been shown that misalignment can affect the shear stress of the wall. It was shown that the gas superficial velocity directly affects the flow behavior, as well as the maximum pressure drop and it standard deviation upstream the sudden expansion of the channel. The experimental investigation using the high-speed photography of gas-liquid flow in a horizontal pipe with a sudden expansion was carried out in [9]. The bubble burst phenomenon is described, and liquid film thickness is estimated.

The PIV and shadowgraph measurements of flow and bubble dynamics in an upward bubbly flow in a square pipe with a sudden expansion were studied in [10]. The experiments were carried out for both laminar and turbulent flow regimes. Authors observed that in a two-phase bubbly flow, smaller bubbles migrate toward the pipe wall and larger ones rise in a tube core region by the effect of steeper velocity gradient with the increasing Reynolds number. The accurate modeling of flow structure, void fraction distribution over a duct cross-section and heat transfer in turbulent bubbly flows with a sudden expansion is very important for safety operations, and the predictions of different emergency situations for the energy equipment of thermal and nuclear power stations.

There are only a few works concerned the numerical study of turbulent bubbly flows behind a pipe or duct with sudden expansion without [11] and with heat transfer [12,13]. The mathematical model for describing a polydispersed bubbly flow in the vertical pipe with sudden expansion and its validation by their own measurements was performed in [8]. The turbulence of the carrier phase was modeled using the *k*-*ω* SST model [14]. The bubble break-up and coalescence dynamics are predicted using the inhomogeneous multiple size group (*h*-MUSIG) model [15].

The mathematical model developed and numerical results of the structure and heat transfer in a polydispersed bubbly turbulent flow downstream of a sudden pipe expansion were presented [12]. The turbulence of the carrier liquid phase is predicted using the model of Reynolds stress transport. Bubble dynamics was predicted taking into account changes in the average volume of bubbles due to the expansion of the change in their density, break-up and coalescence. The study was carried out at the change of initial diameter of air bubbles in the range of *d*<sup>1</sup> = 1–3 mm and their volumetric gas flow rate ratios of *β* = 0%–10%. The effect of the gas volumetric flow rate ratios on the flow structure and heat transfer in the two-phase flow is numerically studied in a bubbly polydispersed flow in a horizontal duct with a single-side BFS [13]. The model based on the Eulerian two-fluid approach.

The modeling of vertical turbulent bubbly flows in a pipe (duct) or in a bubble column requires solving a wide range of time and length scales. The large Eddy simulation (LES) is one of the modern methods for simulating turbulent flows without modeling the fluid phase turbulence. This approach for modeling bubbly flows is much more reasonable than DNS, but it is still too complicated for engineering simulations of turbulent bubbly flows. LES was successfully used for the simulation of two-phase bubbly shear flows [16] or bubble motion in a bubble column [17,18] in the last two decades.

The effect of the bubble fraction and their diameter on the length of the flow recirculating zone, heat transfer and carrier fluid turbulence remain open. The authors did not find works concerning experimental or numerical study of heat transfer in turbulent bubbly flows in BFS. This work is a continuation of our previous papers [12,13] for modeling of flow and heat transfer in the bubbly polydispersed turbulent flows in a pipe [12] or duct [13] with sudden expansion. The main difference of this work from [12] is the method of simulations of the evolution of the spectrum of the bubble diameter. In this paper, the method of *δ*-functions is used [19,20]. Authors [12] used the model of [21] to describe the evolution of bubble size.

Another method for controlling the mean flow patterns, pressure drop, friction and heat transfer rate is the addition of nanoparticles to a backward-facing step flow [22]. The single-phase laminar flow and heat transfer of nanofluid flow in a horizontal backward facing step subjected to a bleeding condition (suction/blowing) is numerically investigated in [22].

The aim of the present paper is the experimental and numerical study of the flow and bubble dynamics, turbulence modification of the carrier fluid phase, mixing and heat transfer enhancement in the upward bubbly flow in a single-side backward-facing step. This paper may be interesting for scientists and research engineers dealing with the problem of flow control and heat transfer enhancement in power equipment.

#### **2. Measurement Setup**

The scheme of the test facility is shown in Figure 1a. The test liquid was stored in the tank (1). A centrifugal pump manufactured by Grundfos (Bjerringbro, Denmark) (2) with a maximum flow rate of 4.2 l/s is used to supply the working fluid. The pump speed and the flow rate of the fluid are set by the PumpMaster PM-P540 frequency converter. To measure the flow rate of the test fluid, an ultrasonic flow meter (3) KARAT-520-32-0 (NPO KARAT, Yekaterinburg, Russia) is used, the measurement uncertainty of the flow rate is 0.2% of the measured value. The test section (duct with sudden expansion) (4) was made of an organic glass with the inner dimensions 1000 × 200 × 20 mm. To create a separation zone, a plexiglass plate is fixed inside the channel. The length and height of the duct before the sudden expansion were L1 = 600 mm and *h*<sup>1</sup> = 8 mm, and after the sudden expansion were L2 = 400 mm and *h*<sup>2</sup> = 20 mm. The step height H = 12 mm, and the expansion ratio was ER = (*H* + *h*1)/*h*<sup>1</sup> = *h*2/*h*<sup>1</sup> = 2.5. A honeycomb was mounted on the duct inlet in front of the plate to establish a uniform fluid flow. The test section is fixed on a frame made of machine-tool aluminum profiles, which provides convenient fastening of the stand elements and the optical system and easy readjustment.

The liquid flow rate through the working section was 0.72 L/s, which corresponds to the mean-mass flow velocity *Um*<sup>1</sup> = 0.55 m/s and Reynolds number based on step height Re*<sup>H</sup>* = *Um*1*H*/ν = 6600. Gas bubbles were introduced into the liquid flow through 9 capillaries with an inner diameter of 0.7 mm located in the lower part of the working section.

The measurements were carried out using the PIV/PLIF system (5), based on the use of a green laser with a constant glow (wavelength 532 nm), full power 1 W and a high-speed camera JET 19 (Kaya Instruments, Haifa, Israel). Since the duct has a large width, measurements in a two-phase flow using the PIV/PLIF method in a situation where the laser knife and the camera's field of view are located perpendicularly is extremely difficult. This is due to the significant overlap of the measurement area by bubbles that move between the region of interest (ROI) and the camera. Therefore, for measurements, the laser knife was wound perpendicular to the wall of the test section, and the camera was directed at an angle of 45 (see Figure 1a). Because of this, we had different scales along the horizontal and vertical axes. The calibration of the linear dimensions was performed for both axes. Measurements were performed along the central axis of the channel with a sudden expansion just behind the BFS (see Figure 1b). The shooting speed during the

experiments was 1000 frames per second. Distilled water with the addition of fluorescent (rhodamine filled) polyamide particles manufactured by Dantec Dynamics (Skovlunde, Denmark) (their diameter was 1–20 μm) was used in measurements. To make bubbles invisible, no Rhodamine was added to the test liquid accorded [23].

**Figure 1.** The scheme of the experimental setup (**a**): 1 is the water tank, 2 is the pump, 3 is the flowmeter, 4 is the test section, 5 is the PIV/PLIF system. (**b**) The scheme of a two-phase bubbly upward flow behind a backward-facing step.

The recognition of bubbles in a separated flow is a nontrivial task and requires the development of new methods and approaches. Examples of images obtained by adding Rhodamine to the flow are shown in Figure 2a–c, obtained at Re = 6600 and b = 0.03 and different sizes of bubbles. The figures show that the bubbles in the separation region of the flow can have a substantially non-spherical shape. The currently developed methods, as a rule, deal with spherical or elliptical bubbles [24]. Development of a deep learning-based image processing technique for bubble pattern recognition and shape reconstruction in dense bubbly flows were performed in [25,26]. In addition, depending on the position relative to the laser knife, the bubbles can be "light" or "dark", as well as cast a shadow and reflect a noticeable glare. In the case without adding a fluorescent dye to the flow, but using fluorescent particles for PLIF, the flow pattern is shown in Figure 2d. An optical threshold filter was used to view the position of particles in the flow. In this case, the recognition of bubbles is not required, however, it is possible to obtain information only on the hydrodynamic characteristics of the liquid phase. However, since the liquid flow structure determines the heat transfer, it may be sufficient for some tasks, including the validation of CFD codes.

A series of 10,000 frames for each regime was accrued. Data processing was carried out by "ActualFlow" software developed in IT SB RAS. Firstly, the subtraction of the mean intensity field averaged over the whole sample range were successively applied to enhance the quality of the raw data. Velocity fields were calculated using the iterative cross-correlation algorithm with a continuous window shift and deformation and 75% overlap of the interrogation windows. The initial size of the interrogation window was chosen to be 64 × 64 pixels but it was subsequently reduced to 32 × 32 pixels. The obtained instantaneous velocity vector fields were then validated with the following procedures: peak validation with the threshold of 2.0 and adaptive median 5 × 5 filter.

**Figure 2.** Single frames of the flow at Re = 6600 and β = 0.03: (**a**) "small" bubbles; (**b**) "medium" bubbles; (**c**) "large" bubbles; (**d**) "small" bubbles without Rhodamine in the test fluid.

In order to validate the proposed method, preliminary measurements of the hydrodynamic structure of the flow were carried out using Laser Doppler anemometry system with a measurement uncertainty ±2%. The deviations of the data obtained using the PIV and LDA in a single-phase flow did not exceed 3%.

#### **3. Mathematical Model and Method of Numerical Realization**

#### *3.1. D RANS + SMC "In-House" Numerical Code*

The flow and heat transfer in the bubbly flow is modeled using the Eulerian approach [26]. It treats the dispersed phase as a continuous medium with properties analogous to those of liquid and this method is widely used for the simulation of turbulent bubbly flow with and without bubble break-up and coalescence processes in various types of flow geometries. The Eulerian approach is based on kinetic equations for a one-point probability density function [27]. The 2D steady-state, incompressible RANS model [12] is used for modeling of turbulent bubbly flow in the BFS. The set of governing equations consists of continuity, two-momentum and energy equations taking into account the effect of bubble presence on the transport processes in the carrier fluid phase. The effect of bubbles on the mean and fluctuational characteristics of turbulent carrier fluid flow is determined by drag, gravity lift, virtual mass, wall lubrication forces, turbulent transport and turbulent diffusion (so called two-way coupling) [12].

The system of 2D RANS equations for the two-phase bubbly flow is given as [23,28]

$$\begin{aligned} \nabla \times (\Phi\_l \rho \mathbf{J} \mathbf{I}) &= 0 \\ \nabla \times (\Phi\_l \rho \mathbf{J} \mathbf{I} \mathbf{I}) &= \Phi\_l (-\nabla P + \rho \mathbf{g}) + \nabla \times (\Phi\_l \rho \mathbf{r}) - \nabla \times \left( \Phi\_l \rho \langle \mathbf{u}' \mathbf{u}' \rangle + \sum\_{k=1}^K \sigma\_k^{R\mathbf{I}} \right) + \\ + (P - P\_{\text{lin}}) \nabla \times \Phi\_l + \sum\_{k=1}^K M\_{\text{lk}} \\ \nabla \times (\Phi\_l \rho \mathbf{C}\_P \mathbf{I} \mathbf{I} \mathbf{T}) &= \nabla \times \Phi\_l \nabla [(\lambda \mathbf{T}) \mathbf{I}] - \nabla \times (\Phi\_l \rho \mathbf{C}\_P \langle \mathbf{u} \theta \rangle) + \sum\_{k=1}^K \frac{\Phi\_{\text{lk}} (\mathbf{T} - \mathbf{T}\_{\text{lk}})}{d\_k} + \\ + \sum\_{k=1}^K \Phi\_{\text{lk}} \mathbf{C}\_{Dk} \mathbf{g}\_{\text{ilt}} \langle \mathbf{u} \theta \rangle \nabla \times \Phi\_l \end{aligned} \tag{1}$$
 
$$\Phi\_l + \Phi\_{\text{b}} = \Phi \equiv 1$$

Hereinafter, the index *k* is partially omitted; *Φ<sup>l</sup>* and *Φ<sup>b</sup>* are the volume fractions of the liquid phase and bubbles, respectively; *K* is the number of groups of bubbles; *P* and *Pin* ≈ *Pb* are the pressures in the liquid phase and on the surface of bubble [1,29], respectively; *σBI* is the influence of the *k*th group of bubbles on the tensor of averaged Reynolds stresses in the liquid phase [1,27]; *τ*, *u u*  and *uθ* are the tensors of viscous stress and Reynolds stress and turbulent heat flow in the carrier phase, respectively; *Ml* = ∑*<sup>K</sup> <sup>k</sup>*=<sup>1</sup> *Mlk* = −*Mb* interfacial term [23,27]; *g* is the gravity acceleration; *gut*,*<sup>k</sup>* is the coefficient of involvement of gas bubbles of the *k*th fraction in the thermal fluctuation movement of the carrier fluid [27]; *τ*Θ*<sup>k</sup>* is the thermal relaxation time. The right-hand side of the momentum equation includes the pressure gradient in carrier phase, gravity, the viscous stress, Reynolds stresses in bubbly flow, interfacial pressure gradient and momentum exchange between carrier fluid and gas phases that arise from the actions from interfacial forces.

The turbulence of the carrier fluid is predicted using the Reynolds stress model (second moment closure—SMC) [30] taking the effect of gas bubbles on fluid phase turbulence [31].

$$\begin{array}{c} \nabla \left( \Phi\_l \rho \boldsymbol{L} \langle \boldsymbol{u}' \boldsymbol{u}' \rangle \right) = \nabla \times \left( \Phi\_l D\_K \right) + \Phi\_l (P\_K + \phi - \varepsilon) + S\_K\\ \nabla \times \left( \Phi\_l \rho \boldsymbol{L} \boldsymbol{\varepsilon} \right) = \nabla \times \left( \Phi\_l D\_\varepsilon \right) + \Phi\_l (P\_\varepsilon - \varepsilon) + S\_\varepsilon\\ \chi - L\_T^2 \nabla^2 \chi = \frac{1}{\varepsilon T\_T} \end{array} \tag{2}$$

Here *<sup>φ</sup>ij* = (<sup>1</sup> − *<sup>k</sup>χ*)*φ<sup>W</sup> ij* + *<sup>k</sup>χφ<sup>H</sup> ij* is the velocity-pressure-gradient correlation, well known as the pressure term [30], *χ* is the blending coefficient, which goes from zero at the wall to 8 unity far from the wall [27]; *φ<sup>H</sup> ij* is "homogeneous" part (valid away from the wall), and *<sup>φ</sup><sup>W</sup> ij* is "inhomogeneous" part (valid in the wall region) and *LT* <sup>=</sup> 0.45 *max k*3/2/*ε*; 80(*ν*/*ε*) 3/4 is the turbulent macro-scale length. Last terms *SK* and *S<sup>ε</sup>* in the right part of the turbulence model take into account the effect of the dispersed phase on transport processes [31].

The dispersed phase is modeled by the Eulerian two-fluid approach. The system of axisymmetric mean transport equations for the *k*th fraction of bubbles has the form [23,28] (to simplify the equations the index *k* is partially omitted).

$$\frac{\frac{\partial \rho\_k}{\partial t} + \nabla \times (a\_{bk} \rho\_{bk} l L\_{bk})}{Dt} = 0$$

$$\frac{D(a\_{bk} \rho\_{bk} l l\_{bk})}{Dt} = a\_{bk} \left(\frac{D\_{bk}}{\tau\_k} \nabla \times (a\_{bk} \rho\_{bk})\right) - \nabla \times (a\_{bk} \rho\_{bk} E \langle u' u' \rangle) + \sum\_{k=1}^{K} M\_{bk} \tag{3}$$

$$\frac{D(a\_{bk} \rho\_{bk} T\_{bk})}{Dt} = -\frac{1}{\tau\_{\rm Qk}} \nabla \times \left(\rho\_{bk} D\_{bk}^{\rm Q} a\_{bk}\right) - a\_{bk} (T - T\_{bk}) \frac{\rho\_{bk}}{\tau\_{\rm Qk}} \rho\_{bk} = a\_{bk} P\_{b} / \left(R\_{b} T\_{bk}\right)$$

The method of *δ*-approximation of [19,20] is used for modeling the polydispersity of the bubbly flow. The break-up and coalescence processes are predicted according to the model of [23,28] and they are described in detail in [32].

The governing transport equations for both gas (1) and dispersed (2) phases, and the turbulence model (3) are solved using a control second-order upwind finite control volumes approach on a staggered grid. The convective terms are discretized using the third order QUICK algorithm and the diffusion terms are numerically solved employing the second-order central difference scheme. The SIMPLEC scheme is used for the coupling of velocity and pressure. The basic grid with 256 × 124 control volumes (CVs) along the longitudinal and transverse directions is used for all numerical simulations. The time derivatives are discretized using the implicit Euler scheme of the first-order accuracy. The components of Reynolds stresses are determined at the same points along the faces of the control volume as the corresponding components of the mean fluid velocity [33]. Grid convergence is verified for three grid sizes is used for the "in-house" code: coarse grid: 128 × 62, basic grid 256 × 124 and fine grid: 400 × 180 control volumes (see Figure 3a).

**Figure 3.** Grid convergence test of 2D RANS + SMC "in-house" code for the bubbly flow at *TW* = 313 K, *T*<sup>1</sup> = *Tb*<sup>1</sup> = 293 K, *β* = 5%, *d*<sup>1</sup> = 3 mm (**a**), single-phase fluid flow 3D RANS+SST (**b**) and LES+WALE (**c**) in a duct with sudden expansion. Re*<sup>H</sup>* = 6600, ER = 2.5.

#### *3.2. D RANS + SST and LES (Ansys CFD Package)*

Three-dimensional RANS and LES predictions for the single-phase fluid flow behind the single-side BFS are performed in the paper using the commercial CFD package Ansys Fluent 2020R2 [34]. Ansys CFD package was used for the additional validation of the "in-house" RANS + SMC numerical code. 3D RANS simulations by Fluent are employed Menter's *k*-*ω* SST model [14] and LES with wall-adapting local eddy-viscosity (WALE) subgrid-scale (SGS) model [35]. The Menter's *k*-*ω* SST model is one of the most popular isotropic two-equation turbulence models for simulation of various turbulent flows. The WALE SGS model is one of a major SGS model for flow simulations in complicated geometries.

The Fluent 3D RANS basic grid consist of 256 × 100 × 100 = 2.56 Mio CVs. Grid convergence is verified for three grid sizes: coarse grid: 150 × 50 × 50 CVs, basic grid: 300 × 150 × 100 and fine grid: 450 × 184 × 150 (see Figure 3b). The control volume method for discretization of the governing equations is used. The computational grid is non-uniform, and more refined grid is performed to all wall surfaces. The grid resolution is chosen optimal after preliminary simulations, and it provides the condition for the position of the first cell from the wall *y*<sup>+</sup> < 1. The profiles of the velocity, turbulent kinetic energy and the rate of their dissipation are set at the entrance to the duct with a backward-facing step. They are simulated in the upstream section with the length 5*H* before the section of flow separation. The wall functions are not used, i.e., the predictions are carried out directly to solid surfaces. Boundary conditions at the duct entrance profiles—velocity, turbulent kinetic energy and their dissipation rate—are set after the preliminary simulations of the flow in the upstream section. The temperature of the incoming water flow is constant *T*<sup>1</sup> = 298 K, and the change in the thermophysical properties is neglected in the predictions.

The LES simulations are conducted using the basic grid with 490 × 124 × 30 = 1.82 Mio CVs. Grid convergence is verified for three grid sizes: coarse grid: 250 × 62 × 20 CVs, basic grid: 490 × 124 × 30 and fine grid: 650 × 190 × 45 (see Figure 3c). It should be noted that the periodic boundary condition in *z* (duct width) direction for the LES method were used. The governing equations are solved numerically using a pressure-based non-iterative time-advancement (NITA-FS) solver with Courant-Friedrichs-Lewy condition CFL < 0.5. The second-order central scheme is used to discretize the convective terms. To predict the flow in the upstream section, the 3D RANS method was conducted. LES modeling is carried out only for a duct with a single-side backward-facing step. The grid resolution is chosen optimal after preliminary simulations, and it provides the condition for the position of the first cell from the wall *y*<sup>+</sup> < 1. The temperature of the incoming water flow was constant *T*<sup>1</sup> = 298 K, and the change in the thermophysical properties was neglected in the calculations. The boundary conditions on the wall behind the step is *TW* = const = 313 K. The end wall of the step is not heated.

Differences in the value of local Nusselt number Nu calculated for the bubbly flow is less than 0.1%. The maximum error *e*max is defined as: *e max i*=1, *N Nu<sup>n</sup> <sup>i</sup>* <sup>−</sup> *Nun*−<sup>1</sup> *i* −6 *max*, where *N* is the total number of CVs in corresponding direction, the subscript *i* is the specific CV and the superscript *n* is the iteration level. The computational grid is nonuniform both in the streamwise and transverse directions. A more refined grid is applied in the recirculation region and in the zones of flow detachment and reattachment and in the inlet region of the duct. The coordinate transformation is suitable for such a two-dimensional problem:

$$
\Delta\psi\_{\bar{j}} = K \times \Delta\psi\_{\bar{j}+1\nu}
$$

where Δ*ψ<sup>j</sup>* and Δ*ψj*-1 are the current and previous steps of the grid in the axial or radial directions and *K* = 1.05 (longitudinal direction) and *K* = 1.03 (transverse direction). The first cell is located at a distance *y*<sup>+</sup> = *yU*\*/*ν* = 0.3–0.5 from the wall surfaces, where *U\** is the friction velocity obtained for the single-phase flow in the inlet of the duct (other parameters being identical). At least 10 CVs were generated to ensure resolution of the mean velocity field and turbulence quantities in the viscosity-affected near-wall region (*y*+ < 10). The time step was Δ*t* = 0.1 ms. The Courant number, which is a necessary condition for the stability of numerical solution of differential equations in partial derivatives, does not exceed 1 for all simulations.

At the inlet streamwise and transverse components, temperatures of the phases and turbulence levels are uniform. No-slip conditions are set on the wall surface for the carrier fluid phase. At the outlet edge, the computational domain condition *∂θ/∂x* = 0 is set for all variables.

#### **4. Comparison with Experimental and Numerical Results**

*4.1. Single-Phase Flow in the BFS*

At the first stage, comparisons were carried out with the measurement data for a single-phase gas duct flow behind a single-side backward-facing step. The well-known experimental results [36] were used for comparisons on the heat transfer. Additionally, comparisons were performed with the semi-empirical correlation data [37] for minimal wall friction coefficient distribution.

The measured [36] and authors' predicted distributions of Stanton number St along the streamwise coordinate are shown in Figure 4a. The measurements were performed in duct with single-side BFS with dimensions: height 188 mm after sudden expansion and width was 500 mm. The step height was *H* = 38 mm and expansion ratio was ER = 1.25. The carrier fluid was air with temperature *T*<sup>1</sup> = 298 K. The mean upstream fluid velocity *Um*<sup>1</sup> = 11 m/s and Reynolds number was Re*<sup>H</sup>* = 2.8 × 104. The Stanton number is calculated using the expression: St = *h*/(*ρCpU <sup>m</sup>*1). Here *h* is the heat transfer coefficient, *ρ* is the density of the carrier fluid flow and *Cp* is the specific heat at constant pressure.

**Figure 4.** Heat transfer coefficient distributions along the streamwise coordinate (**a**) and minimal wall friction values in the reverse zone vs. Reynolds number (**b**). (**a**): Symbols are the experiments of [36], line is our RANS + SMC "in-house" code simulations. (**b**): 1 is the semi-empirical correlation [37], 2 is our RANS+SMC "in-house" code simulations.

The distributions of the minimum value of the wall friction coefficient *Cf,* min in the separation region for single-phase fluid flow vs. the Reynolds number Re*<sup>H</sup>* based on the step height are given in Figure 4b. The measurements of [37] were performed in the duct with single-side backward-facing step with dimensions: length was 1600 mm, height 70 mm after sudden expansion and width was 220 mm. The step height was *H* = 20 mm and expansion ratio was ER = 1.4. The carrier fluid was water with temperature *T*<sup>1</sup> = 298 K. The mean upstream fluid velocity *Um*<sup>1</sup> = 0.024–0.24 m/s and Reynolds number varied Re*<sup>H</sup>* = (0.12—1.22) × 104. The wall shear stress was measured using the electrodiffusion method. Line *1* is the semi-empirical correlation [37] for calculating the wall friction *Cf* ,min = −0.38*ReH*<sup>−</sup>0.57, and points (2) represent the author's predictions. The strong dependency of minimal values of wall friction in the reverse zone on the Reynolds numbers was obtained as well as in experiments [37] and our numerical predictions. Good agreement was obtained between the results [36,37] and the data of authors' numerical predictions.

#### *4.2. Two-Phase Bubbly Flow in Round Vertical Pipe*

Comparison of experimental [38] and authors' predicted data on the distributions of carrier fluid and gas bubbles velocities over the round pipe cross-section is shown in Figure 5a. The mean axial velocities of both carrier fluid (line 2) and gas (line 3) phases in a two-phase flow is always higher than in a single-phase flow (line 1). The velocity profile of gas bubbles (line 3) over the cross section of the pipe has a flatter form in comparison with

the velocity of the fluid (line 2). The mean axial velocity of the bubbles in the upward flow is higher than that one for the liquid phase due to the action of the Archimedes force.

**Figure 5.** Radial profiles of the carrier fluid and gas bubbles (**a**) and mean Sauter gas bubbles diameters (**b**). *β* = 4.6%, Re = *Um*12*R*/*<sup>ν</sup>* = 4.9 <sup>×</sup> 104, *<sup>x</sup>*/(2*R*) = 53.5, 2*<sup>R</sup>* = 50.8 mm, *Um*<sup>1</sup> = 0.986 m/s, *<sup>d</sup>*<sup>1</sup> = 2.4 mm. Symbols are the measurements of [38], curves are the authors predictions. (**a**): 1, 2 are the carrier fluid axial velocities in the single-phase and two-phase bubbly flows respectively; *3* is the gas bubbles axial velocity. (**b**): 1—*Jb* = 0.0473 m/s, *β* = 4.6%, *d*<sup>1</sup> = 2.4 mm; 2—0.113 m/s, 10.3%, 2.5 mm; 3—0.242 m/s, 19.7%, 2.8 mm.

The radial profiles of the mean Sauter bubble diameter along the pipe radius are shown in Figure 5b. An increase in the gas volumetric flow rate ratios causes an increase in the size of bubbles over the pipe cross-section. The bubble size in the near-wall region is slightly higher than the corresponding value for the core part of a two-phase pipe flow. This conclusion is typical both for measurements [36] and for authors' numerical simulations. This due to the fact that small bubbles migrate to the wall due to the action of the lift and turbulent migration forces and they accumulate near the wall. The probability of their coalescence is increased in this case.

#### *4.3. Two-Phase Bubbly Flow in the BFS*

The numerical model was validated against comparison with experimental results for the bubbly flow in a vertical upward pipe with sudden expansion. These results were presented in our previous papers [12,23] and they are not presented here only for brevity. We have carried out the comparison with measurements of [39] and LES of [40] for the isothermal flow regime. The mean-mass velocity was *Um*<sup>1</sup> = 1.78 m/s, which corresponds to the Reynolds number Re*<sup>H</sup>* = 1.1 × 105. The pipe diameter, before the abrupt, was 2*R*<sup>1</sup> = 50 mm, while after expansion, 2*R*<sup>2</sup> = 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 Sauter diameter was *d*<sup>1</sup> = 2 mm. The authors results agree well with experiments and LES data. We have conducted the comparison with present measurements for the mean fluid axial velocity profiles and heat transfer enhancement [23]. The maximal difference between our and other authors measured results and our predictions is up to 20%.

#### *4.4. Two-Phase Bubbly Flow behind a Backward-Facing Step. Our Measurements and Numerical Simulation and Their Discussion*

Based on a successful comparison of the results obtained using different methods, studies of the local structure of the two-phase bubbly flow in a BFS are carried out. Symbols and curves are the authors' measurements and predictions, respectively. The transverse profiles of measured and predicted mean streamwise fluid velocities at the Re*<sup>H</sup>* = 6600 in the single-phase fluid flow (lines 1) and at gas volumetric flow rate ratio *β* = 3% (lines 2) along the duct length are shown in Figure 6. The carrier fluid velocity was normalized to

the maximum value of liquid velocity. The point of flow reattachment in the single-phase fluid flow is located at a distance of (7–8)*H*, which is in good agreement with the literature data [41,42]. The following characteristic features of the bubbly flow can be distinguished. Immediately behind the step, the flow recirculation zone is observed. The gas bubbles badly penetrate this region and the effect of gas bubbles on the flow is rather weak in this area. The position of the flow reattachment point shifts towards the step for the two-phase flow and, in general, the establishment of a profile in the channel after the expansion occurs faster in a two-phase flow than in a single-phase flow, which was previously reported, for example, in [23]. Our numerical results capture the main features of the bubbly flow behind the BFS described above. The main difference between the experimental and numerical results for the two-phase bubbly is obtained in the near-wall zone in the flow relaxation region.

**Figure 6.** The transverse profiles of the carrier phase mean streamwise velocity along the duct length. Symbols and curves are the authors' measurements and predictions respectively. Re*<sup>H</sup>* = 6600, ER = 2.5, *β* = 3%. 1—single-phase flow (*β* = 0), 2—*β* = 3%. (**a**)—*x/H =* 2, (**b**)—4, (**c**)—6, (**d**)—8, (**e**)—10.

#### **5. Numerical Results and Its Discussion**

#### *5.1. Single-Phase Fluid Flow in a Backward-Facing Step Flow*

The distribution of mean streamwise fluid velocity component at various distances from the BFC is presented in Figure 7. The measurements are conducted using the PIV. The 3D RANS + SST predicted contours of the longitudinal velocity component and streamlines at a few stations along the *z* axis are shown in Figure 8a. The computations are carried out using the commercial CFD package Ansys Fluent 2020R2. This figure clearly shows the three-dimensional flow pattern in the recirculation area. The length of the recirculation zone is maximum in the central section (*z*/*H* = 8.3). The length of the recirculation zone behind the backward-facing step increases with distance from the side walls of the duct. When flowing around the backward-facing step, the flow detaches the sharp edge, and a separation bubble is formed. A long recirculation zone and, characteristic of it, a secondary vortex is formed behind the edge. It has the opposite direction of flow rotation to the main one, and the rotation velocity in it is relatively low, which leads to the formation of a stagnant zone in this region.

**Figure 8.** The 3D RANS + SST predicted contours of the longitudinal velocity component and streamlines (**a**) and turbulent kinetic energy (**b**) at the *z* = 10, 50, 100, 190 mm (*z*/*H* = 0.8, 4.2, 8.3 and 15.8).

The change in the turbulent kinetic energy behind the step in various sections along the *z* axis is presented in Figure 8b. A turbulization of the flow occurs when flowing in the BFS due to additional vortex formation. The largest values of TKE are located in the region of the mixing shear layer, and their minimum values are observed in the region of the secondary vortex. The turbulent kinetic energy increases toward the center duct. The largest values of turbulent kinetic energy in the region are observed in the central section.

The transverse profiles of single-phase flow mean streamwise velocity behind a singleside backward-facing step are shown in Figure 9. Here symbols and curves are the authors' measurements and predictions, respectively. The first three sections are located in the recirculation region, the fourth station is set close to the reattachment point and last two sections are set in the flow relaxation zone. It was obtained a qualitatively agreement between authors' measurements and LES and RANS results.

**Figure 9.** The transverse profiles of single-phase mean longitudinal velocity behind a backwardfacing step. Symbols and curves are the authors' measurements and predictions, respectively. Re*<sup>H</sup>* = 6600, ER = 2.5.

#### *5.2. Upward Bubbly Flow in a Backward-Facing Step Flow*

Numerical simulations are carried out for a fluid flow of a water and polydispersed air bubbles at atmospheric pressure and all these results given in this subsection are obtained using the "in-house" 2D RANS + SMC model. The upward bubbly flow in the duct with a single-side backward-facing step is considered (see Figure 1b). The Reynolds number Re = *Um*1*H*/ν = 6600, and the mean-mass velocity at the inlet *Um*<sup>1</sup> = 0.55 m/s. The height of the duct before separation *h*<sup>1</sup> = 8 mm, height of the duct after separation *h*<sup>2</sup> = 20 mm, and the step height *H* = 12 mm, the expansion ratio ER = (*H* + *h*1)/*h*<sup>1</sup> = *h*2/*h*<sup>1</sup> = 2.5, the wall temperature *TW* = const = 313 K and the initial temperatures of carrier liquid and gas bubbles *T*<sup>1</sup> = *Tb*<sup>1</sup> = 293 K. In the inlet section, the mean phase velocities have the same value. The initial gas volumetric flow rate ratios are varied *β* = 0–10% and the mean initial Sauter diameter of air bubbles *d*<sup>1</sup> = 3 mm. All simulations are performed for four monodispersed *δ*-functions (modes) of air bubbles at the inlet cross-section: 1—*d*/*d*<sup>1</sup> = 0.33, 2—0.66, 3—1, and 4—1.33. The volume fraction of each fraction is *Φ*1(1 mm) = 0.031Φ, *Φ*2(2 mm) = 0.12Φ, *Φ*3(3 mm) = 0.6Φ, and *Φ*4(4 mm) = 0.15Φ, where *Φ* is the total volume fraction of gas bubbles. Authors did not carry out the measurements of the distributions of the diameters of gas bubbles at the inlet cross-section. Therefore, the choice of just such values of gas bubbles diameters and their volume fractions is based on our preliminary numerical predictions. There is no steam generation on wall surface. The same assumptions were used in our previous recent numerical simulations [12,13,29], but they may be important in other thermal boundary conditions on the wall [43].

In Figure 10 the transverse profiles of the total local void fraction (a) and total bubble diameter (b) of all monodispersed four modes are shown along the duct length in the upward bubbly flow behind the backward-facing step. The so-called "top-hat" distribution [44] of the local void fraction is obtained at the inlet (line 1) (the cross-section of sudden expansion). This is one of the most typical distributions of void fraction in a vertical upward bubbly flow for a low gas bubble concentration. The obvious maximum of void fraction in the near-wall regions is typical for this type of void fraction distribution [44]. Further, there is a significant change in the predicted distributions of the local void fraction over the duct cross-section. It can be explained by a significant increase in the two-phase flow cross-section behind the BFS (see Figure 10a). The "right" local maximum of the void fraction is also located near the "right" wall of the duct, but its value is noticeably less than for the inlet section (lines 2–5). The "left" local maximum of the void fraction

becomes much less noticeable along the duct length by bubbles scattering and diffusion across the duct cross-section. In the area of the recirculation zone of the two-phase flow (curves 2–4), an extremely small number of gas bubbles are obtained. It can be explained by their relatively large initial diameter. It was recently shown in our previous work [12] that large bubbles practically do not penetrate the recirculation region in a pipe with sudden expansion and are observed mainly in the mixing layer and in the flow core. The bubbles distribution after the reattachment point is characterized by the presence of bubbles over the entire all cross section of the duct by the effect forces acting on a bubble in a transverse direction and low maximum in the near-wall zone of the duct (line 5). The profile of local void fraction at *x*/*H* = 10 has almost uniform distribution. The profiles of the local void fraction are characterized by a zero-value close to the duct wall. This is explained by the fact under the action of the wall force a bubble cannot approach a wall. The profiles of dimensionless bubble diameter of all monodispersed four δ-functions are presented in Figure 10b. The maximum bubbles size is obtained in the near-wall zone of the duct. The bubbles diameter in the recirculation region is less than that one in the flow core.

**Figure 10.** The transverse profiles of total local void fraction (**a**) and bubble diameter (**b**) along the duct length in the upward bubbly flow behind the backward-facing step. Re*<sup>H</sup>* = 6600, ER = 2.5, *β* = 3%. (**a**): 1—*x*/*H* = 0, 2—2, 3—4, 4—6, 5—10; (**b**): 1—*x*/*H* = 2, 2—4, 3—6, 4—10.

The "in-house" 2D RANS + SMC numerical results for the sum (total) of all four fractions of air bubbles were shown in Figure 11. To better understand the complicated and coupled break-up and coalescence processes in bubbly turbulent flows, the results of local void fractions and bubble diameter distributions obtained for four various bubble fractions are important, and these results are presented in Figure 11. The cross-section *x*/*H* = 4 is located in the area of flow recirculation. The void fraction profiles of the smallest (line 1) and middle (curve 2) bubbles have only one local maximum located close to the duct wall. Only the bubbles of the smallest size (line 1) are observed in the recirculation zone of the duct (see Figure 11a). The void fraction profiles of large bubbles (curves 3 and *4*) have two local maxima and are set close to the duct wall and in the shear mixing layer. Two obvious local maxima are revealed in the distribution of void fraction of all modes (line 5) too.

**Figure 11.** The transverse profiles of local void fractions (**a**) and bubble diameter (**b**) for various bubble size modes. Re*<sup>H</sup>* = 6600, ER = 2.5, *β* = 3%, *x*/*H* = 4. (**a**): 1—*d*/*d1* = 0–0.33, 2—0.33–0.66, 3—0.66–1, 4—>1, 5—sum of all modes; (**b**): 1—*d* = 0–0.33, 2—0.33–0.66, 3—0.66–1, 4—>1.

The distribution of all fractions of bubble diameters across the duct cross-section is practically uniform (see Figure 11b). Bubbles of the smallest investigated diameter *d*/*d1* ≤ 0.33 (line 1) are found over the entire section of the channel, both in the flow recirculation zone and in the flow core. The smallest air bubbles can come closer to the duct wall than the larger bubbles due to the effect of transverse forces (lift, turbulent migration, turbulent diffusion and wall force). The largest bubbles *d*/*d1* ≥ 1 (line 5) are located mainly in the near-wall part of the channel, which confirms the data of our numerical predictions presented in Figure 11a,b. Bubbles of the other two monodispersed fractions (*d*/*d1* = 0.33–0.66, 0.66–1, lines 2 and 3) also practically do not penetrate the flow separation region.

The effect of the volumetric flow rate gas content on the profiles of the longitudinal liquid velocity is shown in Figure 12. The carrier fluid (water) velocity in the two-phase bubbly flow (curves 2 and 3) is slightly higher that one in the single-phase flow (line 1). This effect increases with an increase in the concentration of air bubbles.

**Figure 12.** Profiles of mean axial fluid velocity in bubbly flow in the backward-facing step. 1—*β* = 0 (single-phase fluid flow, other conditions being identical), 2—2%, 3—5%.

The effect of the gas volumetric flow rate ratios (a) and gas bubbles diameter at the inlet (b) on the carrier fluid phase turbulent kinetic energy is shown in Figure 13. Carrier fluid phase turbulence level for a 2D flow is estimated using the expression: <sup>2</sup>*<sup>k</sup>* = *uiui* = *<sup>u</sup>*<sup>2</sup> + *<sup>v</sup>*<sup>2</sup> + *<sup>w</sup>*<sup>2</sup> ≈ 1.5 *u*<sup>2</sup> + *v*<sup>2</sup> . The maximum of TKE of the fluid in the two-phase flow is located in the shear mixing layer. The significant enhancement of the carrier fluid phase turbulence level in the two-phase flow (up to 40% at *β* = 10%) is shown in comparison with the single-phase flow (bold solid curves). The increase of initial bubble size causes a significant increase in the turbulent kinetic energy in the bubbly flow (up to 25% at *d*<sup>1</sup> = 3 mm). The additional production of carrier fluid phase turbulence is explained by vortex formation upon streamlining of the gas bubbles by the carrier fluid flow. The profiles of TKE at a small value of *β* agree qualitatively with those for the case of one-phase flow in a duct with single-side BFS. The same tendencies were obtained in our previous works [12,23] for a pipe with sudden expansion.

**Figure 13.** The effect of the gas volumetric flow rate ratios (**a**) and gas bubbles diameter at the inlet (**b**) on the carrier fluid phase turbulence.

In Figure 14 the effect of the addition of the gas bubbles on heat transfer from the wall to the two-phase bubbly flow. The local convective Nusselt number at constant wall temperature is estimated:

$$\text{Nu} = -(\partial T / \partial y)\_W H / (T\_W - T\_m)\_s$$

where (*∂T/∂y*) *<sup>W</sup>* is the gradient of the fluid temperature on the wall, *TW* is the wall temperature and *Tm* is the mean temperature of the carrier fluid (water) in this section. Lines 1 are Nusslet number in the fully developed single-phase duct flow and dashed

curves 2 are the predictions of heat transfer in the single-phase duct flow with sudden expansion for the fluid (water) flow with other conditions being identical.

**Figure 14.** Heat transfer in the bubbly flow behind the backward-facing step along the streamwise coordinate. (**a**): *d*<sup>1</sup> = 3 mm, 1—Nusslet number value in the fully developed single-phase duct flow, 2—*β* = 0 (single-phase fluid flow, other conditions being identical), 3—2%, 4—5%; (**b**): *β* = 5%, 1—Nusslet number value in the fully developed single-phase duct flow, 2—*d*<sup>1</sup> = 0 (single-phase fluid flow, other conditions being identical), 3—0.5 mm, 4—2 mm, 5—3 mm.

A noticeable increase of heat transfer (up to 35% at *d*<sup>1</sup> = 3 mm in comparison with the single-phase fluid flow) with an increase in the volumetric flow rate ratio *β* is observed (see Figure 14a). This is explained by an increase in the velocity and temperature gradients and turbulization of the carrier fluid (water) phase in the near-wall region of the duct. The position of the peak of heat transfer shifts upstream and at *β* = 5% is *x*Nu\_max/*H* = 4.5, and the length of the recirculation zone approximately coincides with it, *xR*/*H* = 4.8 at *β* = 5%. The same values for the single-phase fluid flow are *x*Nu\_max/*H* = 6.8 and *xR*/*H* = 7. The position of the peak of heat transfer rate *x*Nu\_max is close to the position of the flow reattachment point *xR* in two-phase bubbly flow. The same tendency was observed in [41,42] for single-phase fluid flows in a BFS. The effect of bubble diameter in the inlet on the Nusselt number distributions along the duct length is given in Figure 14b. The largest increase in heat transfer (approximately 35% at fixed value *β* = 5%) is characteristic for the largest gas bubbles with the initial diameter *d*<sup>1</sup> = 3 mm.

The effect of gas volumetric flow rate ratios *β* on values of minimal wall friction coefficient *Cf*,min/*Cf*,min,0, maximal heat transfer enhancement ratios Numax/Numax,0, recirculation length *xR*, maximal values of turbulence modification ratios *k*max and position of heat transfer maximum *x*max are shown in Figure 15. All these variables are normalized on the value in the single-phase fluid flow and subscript "0" is the in the single-phase flow parameter with other conditions are identical to the two-phase bubbly flow.

It is obtained that the magnitude of minimal wall friction coefficient decreases with the increase in the bubble concentration and for *β* = 10% the minimal value of wall friction ratio is *Cf*,min/*Cf*,min,0 ≈ 0.75. The reduction of minimal value of wall friction ratio is small (up to 10%) for the *β* ≤ 5%. Almost linear convective heat transfer enhancement is observed for the range studied of gas volumetric flow rate ratios and for *β* = 10% the heat transfer augmentation ratio is Numax/Numax,0 ≈ 1.75. The addition of gas bubbles leads to the significant shortening of the recirculation length ratio of bubbly flow. The shortening of the recirculating zone is almost twice in comparison with the case of the single-phase fluid flow. The main reason is the flow turbulization by flowing around the gas bubbles. The maximal value of turbulent kinetic energy modification is up to 50% at *β* = 10%. It causes the intensification of mixing process. The same trends were obtained for the single-phase separated flows [41,42]. It is well known [41,42] that the position of maximal value of heat transfer is close to the flow reattachment point for the single-phase for both gas (air) and fluid (liquid) flows behind a BFS and a sudden pipe expansion. The addition of relatively large air bubbles shifted the position of maximal value of heat transfer far from the reattachment point for the bubbly flow in the pipe with sudden expansion [12]. In the case of bubbly flow in the backwards-facing step we observe the same trends. The position of the heat transfer maximum is located after the reattachment point.

**Figure 15.** The magnitudes of minimal wall friction coefficient (1), maximal heat transfer enhancement ratios (2), recirculation length (3), maximal values of turbulence modification ratios (4) and position of heat transfer maximum (5) depending on gas volumetric flow rate ratios *β*. Re*<sup>H</sup>* = 6600, ER = 2.5.

#### **6. Conclusions**

The PIV/PLIF measurements of mean streamwise flow structure and RANS numerical simulation of the bubbly polydispersed upward duct flow in a backward-facing step are performed. Measurements of the carrier fluid phase velocity are performed using the PIV/PLIF system. The numerical model is based on the Eulerian two-fluid approach. Turbulence of the carrier fluid phase is predicted using the Reynolds stress model.

Experimental study of the characteristics of the motion of gas bubbles behind a single side backward facing step with variations in the volumetric flow rates of carrier fluid and gas have been carried out. The distribution of the mean streamwise gas bubbles and liquid velocities are obtained. It is shown that the bubbles slow down at a distance from the flow detachment cross-section, and it is associated with the slowing down of the liquid velocity due to the increase in the duct cross-section. The bubbles have a more curvilinear trajectory in the separation zone than in the duct before the sudden expansion cross-section. It is shown that bubble clusters can form in the flow relaxation region.

Small bubbles are presented over the entire duct cross-section, and the larger bubbles mainly observed in the shear mixing layer and flow core. The recirculation length in two-phase bubbly flow is shorter (up to twice) than that one in the case single-phase flow. The position of the heat transfer maximum is located after the reattachment point. The effect of the gas volumetric flow rate ratios on the flow and heat transfer in the two-phase flow is numerically studied. The addition of air bubbles results in a significant increase in the heat transfer (up to 75%).

**Author Contributions:** Conceptualization, P.D.L. and M.A.P.; methodology, T.V.B., A.V.C., P.D.L. and M.A.P.; investigation, T.V.B., I.A.E., D.V.K., P.D.L. and M.A.P.; data curation, M.A.P.; formal analysis, P.D.L. and M.A.P.; writing—original draft preparation, P.D.L. and M.A.P.; writing—review and editing, P.D.L. and M.A.P.; resources, P.D.L. and M.A.P.; project administration, P.D.L. and M.A.P.; All authors have read and agreed to the published version of the manuscript.

**Funding:** The numerical part of the paper was funded by Ministry of Science and Higher Education of the Russian Federation (mega-grant 075-15-2021-575) and the measurements were performed under IT SB RAS state assignment (Project No. 121032200034-4).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are contained within the article.

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


#### **References**

