RANS Modeling of Turbulent Flow and Heat Transfer in a Droplet-Laden Mist Flow through a Ribbed Duct

: The local structure, turbulence, and heat transfer in a ﬂat ribbed duct during the evaporation of water droplets in a gas ﬂow were studied numerically using the Eulerian approach. The structure of a turbulent two-phase ﬂow underwent signiﬁcant changes in comparison with a two-phase ﬂow in a ﬂat duct without ribs. The maximum value of gas-phase turbulence was obtained in the region of the downstream rib, and it was almost twice as high as the value of the kinetic energy of the turbulence between the ribs. Finely dispersed droplets with small Stokes numbers penetrated well into the region of ﬂow separation and were observed over the duct cross section; they could leave the region between the ribs due to their low inertia. Large inertial droplets with large Stokes numbers were present only in the mixing layer and the ﬂow core, and they accumulated close to the duct ribbed wall in the ﬂow towards the downstream rib. An addition of evaporating water droplets caused a signiﬁcant enhancement in the heat transfer (up to 2.5 times) in comparison with a single-phase ﬂow in a ribbed channel.


Introduction
The intensification of heat transfer in the internal cooling channels of gas turbine (GT) blades remains one of the key problems due to the constant growth in the inlet gas temperature of the GT. This temperature already reaches 2000 K and significantly exceeds the allowable temperatures for the long-term operation of the blades and power equipment of gas turbines [1][2][3][4]. Therefore, cooling the working surfaces of heat-loaded elements is an important and urgent problem of heat transfer. Various cooling methods (film cooling, jet impingement cooling, internal convective cooling, thermal barrier coatings, and spray cooling by the evaporation of various atomized droplets) have been developed for the effective thermal protection of working surfaces and increasing the operating times of power equipment elements. Internal convective cooling is a reliable and simple method for efficient cooling and heat removal from the GT heat-loaded elements.
One of the most effective methods for increasing heat transfer is the use of passive heat transfer intensifiers with various surface shapes. The use of various ribs or obstacles installed on a duct wall is one of the most effective ways to increase heat transfer (see monographs [5][6][7]). The rib height, h; duct height, H; rib pitch, p; obstacle shape; rib-tochannel height expansion ratio, ER = h/H; pitch-to-height ratio, p/h; and some other factors have a great effect on the formation and development of the recirculation region and heat transfer in such flows.
The heat transfer enhancement (HTE) of ribbed ducts (by 2-5 times) is accompanied by a significant increase in the pressure drop (of more than ten times) for most of these surfaces [1,3,4]. Two-dimensional obstacles most often have the form of ribs and protrusions of various configurations located at different angles to the flow on the duct walls [4][5][6][7]. They deflect and mix the flow, give rise to multiscale separated flows, and generate vorticity

Mathematical Models
The paper considers the flow dynamics and heat transfer in 2D two-phase gas-droplet turbulent flow in the presence of interfacial heat transfer between the ribs. The two-fluid Euler approach is used to describe the flow dynamics and heat and mass transfer in the gaseous and dispersed phases [26,27]. The carrier phase turbulence is predicted using the elliptical Reynolds stress model [28], taking into account the effect of droplets [29,30]. The dispersed phase (water droplets) is described using steady-state continuity equations, two momentum equations, and energy equations. The authors used their own in-house code for all numerical simulations presented in this paper.

Governing Equations for the Two-Phase Turbulent Mist Phase
The set of incompressible steady-state 2D RANS equations of the carrier phase includes continuity equations, two momentum equations (in streamwise and transverse directions), energy equations, and steam diffusion into the binary air-steam medium [25]. The effect of evaporating water droplets on the motion and heat transfer in the carrier phase (air) is considered using the sink or source terms.
Here, U i (U x ≡ U, U y ≡ V) and u i (u x ≡ u , u y ≡ v ) are components of mean gas velocities and their pulsations; x i are projections on the coordinate axis; 2k = u i u i = u 2 + v 2 + w 2 ≈ u 2 + v 2 + 0.5 u 2 + v 2 ≈ 1.5 u 2 + v 2 is the kinetic energy of gasphase turbulence; τ = ρ L d 2 /(18ρνW); W = 1 + Re 2/3 L /6 is the particle relaxation time, taking into account the deviation from the Stokes power law; and Re L = |U − U L |d/ν is the Reynolds number of the dispersed phase.
The turbulent heat u j t = − ν T Pr T ∂T ∂x j , and the mass u j k V = − ν T Sc T ∂K V ∂x j fluxes in the gas phase are predicted using simple eddy diffusivity (Boussinesq hypothesis). The constant value of the turbulent Prandtl and Schmidt numbers, Pr T and Sc T equal to 0.9, is used in this work.

Evaporation Model
The set of Eulerian Equation (1) of two-phase flow is supplemented by the equations of heat transfer on the droplet surface and the conservation equation of steam on the surface of the evaporating droplet [31]. It is assumed that the temperature over the droplet radius is constant [31].
Here, λ L is the coefficient of heat conductivity of the droplet; α and α P are the heat transfer coefficient for the evaporating droplet and non-evaporating particle, respectively; T L is the temperature of the droplet; J is the mass flux of steam from the surface of the evaporating droplet; L is the latent heat of evaporation; ρ is the density of the gas-steam mixture; D is the diffusion coefficient; and K * V is the steam mass fraction at the "steamgas mixture-droplet" interface, corresponding to the saturation parameters at droplet temperature T L . Subscript "L" corresponds to the parameter on the droplet surface. The Jacob number, Ja = C P (T − T L )/L, is the ratio of sensible heat to latent heat during droplet evaporation. It characterizes the rate of the evaporation process and is the reciprocal of the Kutateladze number, Ku. For our conditions, the Jakob number is Ja ≤ 0.01.
The expression for the diffusional Stanton number has the form We can insert Equation (4) into Equation (3). Equation (3) can be written in the form is the diffusion parameter of vapor (steam) blowing, determined with the use of a saturation curve.
A droplet evaporates at the saturation temperature, and the temperature distribution inside a droplet is uniform. The droplet temperature along the droplet radius remains constant because the Biot number is Bi = α L d 1 /λ L << 1 and the Fourier number is Fo = τ eq /τ evap << 1. Here, τ eq is the period when an internal temperature gradient inside a droplet exists, and τ evap is the droplet's lifetime. In this case, a droplet evaporates at the saturation temperature, and the temperature distribution inside a droplet is uniform.

The Elliptic Blending Reynolds Stress Model (RSM) for the Gas Phase
In the present study, the low-Reynolds-number elliptic blending RSM of [28] is employed. The transport equations for u i u j and the kinetic energy dissipation rate, ε, are written in the following general form: Here, P ij is the stress-production term, T T and L T are the turbulent time and geometrical macroscales, and φ ij is the velocity-pressure-gradient correlation, well-known as the pressure term. The blending model (8) presented in [28] is used to predict φ ij in Equations (6) and (7), where β is the blending coefficient, which goes from zero at the wall to unity far from the wall; φ H ij is the "homogeneous" part (valid away from the wall) of the model; and φ W ij is the "inhomogeneous" part (valid in the wall region). The other constants and functions of the turbulence model are presented in detail in [28]. The last terms of the system of Equations (6) and (7), A L and ε L , represent the effects of particles on carrier phase turbulence [29,30].

Governing Equations for the Dispersed Phase
The set of incompressible steady-state 2D governing mean equations for the dispersed phase consists of continuity equation, two momentum equations (in streamwise and transverse directions), energy equations.
Water 2022, 14, 3829 where D Lij and D Θ Lij are the turbulent diffusivity tensor and the particle turbulent heat transport tensor [29,30], τ Θ = C PL ρ L d 2 /(12λY) is the thermal relaxation time, and Y = 1 + 0.3Re 1/2 L Pr 1/3 . The set of governing mean equations for the dispersed phase (8)(9)(10) is completed by the kinetic stress equations, temperature fluctuations, and turbulent heat flux in the dispersed phase, which are in the form presented in [29,30].
The volume fraction of the dispersed phase is lower (Φ 1 < 10 −4 ), and the droplets are finely dispersed (d 1 < 100 µm); therefore, the effects of interparticle collisions and break-up are neglected [25,32,33]. Droplet bag break-up is observed at We = ρ(U S − U L ) 2 d/σ ≥ We cr = 7 [33]. Here, U S = U + u S and U L are the gas velocity seen by the droplet [34] and the mean droplet velocity, respectively, U is the mean gas velocity (derived directly from the RANS predictions), u S is the drift velocity between the fluid and the particles [34], and ρ and ρ L are the densities of the gas and dispersed phases. For all droplet sizes investigated in the present paper, the Weber number is very small (We << 1). Droplet fragmentation at its contact with a duct wall also is not considered. The effect of break-up and coalescence in the two-phase mist flow can be neglected due to a low droplet volume fraction at the inlet (Φ 1 = M L1 ρ/ρ L < 2 × 10 −4 ). Here, M L1 is the initial droplet mass fraction, and ρ L is the density of the dispersed phase.
A scheme of the flow is shown in Figure 1. A similar Euler approach was used by the authors to describe gas-droplet axisymmetric flows behind a sudden pipe expansion [25] and behind a backward-facing step in a flat duct [35]. transport tensor [29,30], ( ) The set of governing mean equations for the dispersed phas the kinetic stress equations, temperature fluctuations, and turbu persed phase, which are in the form presented in [29,30].
The volume fraction of the dispersed phase is lower (Φ1 < 10 finely dispersed (d1 < 100 μm); therefore, the effects of interparti up are neglected [25,32,33]. Droplet bag break-up is observed at W = 7 [33]. Here, = + ' S S U U u and UL are the gas velocity seen by mean droplet velocity, respectively, U is the mean gas velocity (d

RANS predictions),
' S u is the drift velocity between the fluid a ρ and ρL are the densities of the gas and dispersed phases. For all d in the present paper, the Weber number is very small (We << 1). D its contact with a duct wall also is not considered. The effect of b in the two-phase mist flow can be neglected due to a low drople inlet (Φ1 = ML1ρ/ρL < 2 × 10 −4 ). Here, ML1 is the initial droplet ma density of the dispersed phase. A scheme of the flow is shown in Figure 1. A similar Euler ap authors to describe gas-droplet axisymmetric flows behind a sud and behind a backward-facing step in a flat duct [35].

Numerical Solution
The solution was obtained using the finite volume method on staggered grids. The QUICK procedure of the third order of accuracy was used for the convective terms. Central differences of the second order of accuracy were used for diffusion fluxes. The pressure field was corrected according to the agreed finite volume SIMPLEC procedure. The components of the Reynolds stress of the carrier fluid phase were simulated according to the method proposed in [36]. The components of the Reynolds stress were determined at the same points along the control volume faces as the corresponding components of the average velocity of the carrier phase. The computational grid consisted of rectangular cells. It was inhomogeneous and thickened towards all solid walls, which was necessary to resolve the details of the turbulent flow in the near-wall zone (see Figure 2). In the viscous sublayer, at least 10 computational volumes (CVs) were set. The correct simulation of sharp gradients of two-phase flow parameters was necessary. The coordinate transformation given in [37] was suitable for such a two-dimensional boundary layer problem. the same points along the control volume faces as the corresponding components of the average velocity of the carrier phase. The computational grid consisted of rectangular cells. It was inhomogeneous and thickened towards all solid walls, which was necessary to resolve the details of the turbulent flow in the near-wall zone (see Figure 2). In the viscous sublayer, at least 10 computational volumes (CVs) were set. The correct simulation of sharp gradients of two-phase flow parameters was necessary. The coordinate transformation given in [37] was suitable for such a two-dimensional boundary layer problem. All predictions were carried out on a "medium" grid containing 256 × 120 control volumes (CVs). The first computational cell was located at a distance from the wall of y+ = u*y/ν ≈ 0.5 (the friction velocity u* was determined for a single-phase air flow with other identical parameters). Additionally, simulations were carried out on grids containing "coarse" 128 × 60 and "fine" 512 × 200 CVs. The difference in the results of the calculations of the wall friction coefficient (a) and the Nusselt number (b) for two-phase flow did not exceed 0.1% (see Figure 3). The Nusselt number at TW = const was determined by the formula: where TW and Tm are the wall and the mass-averaged temperatures of the gas in the corresponding cross section. All predictions were carried out on a "medium" grid containing 256 × 120 control volumes (CVs). The first computational cell was located at a distance from the wall of y + = u * y/ν ≈ 0.5 (the friction velocity u * was determined for a single-phase air flow with other identical parameters). Additionally, simulations were carried out on grids containing "coarse" 128 × 60 and "fine" 512 × 200 CVs. The difference in the results of the calculations of the wall friction coefficient (a) and the Nusselt number (b) for two-phase flow did not exceed 0.1% (see Figure 3). The Nusselt number at T W = const was determined by the formula: where T W and T m are the wall and the mass-averaged temperatures of the gas in the corresponding cross section.  Periodic boundary conditions were set at the inlet of the computational domain. Initially, a single-phase fully hydrodynamically developed air flow was supplied to the inlet to the computational domain L0 = 10p, where p is the rib pitch (the spacing between upstream and downstream ribs). The 1st rib was installed at the end of this domain. The output parameters from section L0 were the input values for section L1 = 10p, located between the 1st and 2nd ribs (see Figure 1). All simulations were performed for the twodimensional case of a gas-droplet flow for the 2nd and 3rd obstacles. Drops were fed into a single-phase turbulent air flow along the entire cross section of the duct in the inlet cross section behind the 2nd rib. The initial temperatures of the gas and dispersed phases at the inlet to the computational domain were T1 = TL1 = 293 K. The boundary condition TW = const = 373 K was set on the ribbed wall; the opposite smooth (without obstacles) wall of the flat duct was adiabatic. The entire ribbed duct surface and all the ribs were heated to eliminate the influence of the possible formation of liquid spots during the deposition of droplets on the wall from a two-phase mist flow. The impermeability and no-slip conditions for the gas phase were imposed on the duct walls. For the dispersed phase on the duct wall, the boundary condition of the "absorbing wall" [30] was used when a droplet did not return to the flow after contact with the wall surface. All droplets deposited from two-phase flow onto the wall momentarily evaporated. Thus, the pipe surface was always dry, and there was no liquid film or spots of deposited droplets formed on the wall Periodic boundary conditions were set at the inlet of the computational domain. Initially, a single-phase fully hydrodynamically developed air flow was supplied to the inlet to the computational domain L 0 = 10p, where p is the rib pitch (the spacing between upstream and downstream ribs). The 1st rib was installed at the end of this domain. The output parameters from section L 0 were the input values for section L 1 = 10p, located between the 1st and 2nd ribs (see Figure 1). All simulations were performed for the twodimensional case of a gas-droplet flow for the 2nd and 3rd obstacles. Drops were fed into a single-phase turbulent air flow along the entire cross section of the duct in the inlet cross section behind the 2nd rib. The initial temperatures of the gas and dispersed phases at the inlet to the computational domain were T 1 = T L1 = 293 K. The boundary condition T W = const = 373 K was set on the ribbed wall; the opposite smooth (without obstacles) wall of the flat duct was adiabatic. The entire ribbed duct surface and all the ribs were heated to eliminate the influence of the possible formation of liquid spots during the deposition of droplets on the wall from a two-phase mist flow. The impermeability and no-slip conditions for the gas phase were imposed on the duct walls. For the dispersed phase on the duct wall, the boundary condition of the "absorbing wall" [30] was used when a droplet did not return to the flow after contact with the wall surface. All droplets deposited from two-phase flow onto the wall momentarily evaporated. Thus, the pipe surface was always dry, and there was no liquid film or spots of deposited droplets formed on the wall [25,31,35]. This assumption for the heated surface is valid (see, for example, papers [25,35]). Furthermore, this condition is valid if the temperature difference between the wall and the droplet is greater than T W − T L ≥ 40 K [38]. In the outlet cross section, the conditions for the equality to zero of the derivatives of all variables in the streamwise direction were set.

Model Validation
At the first stage, a comparison with the data of recent LDA measurements [39] for a single-phase air flow in the presence of ribs was performed. The results of the experiments and predictions are shown in Figure 4. This figure shows comparisons of measured and predicted data in the form of transverse profiles of mean longitudinal velocity, U/U m1 (a), and the velocity of its fluctuations, u'/U m1 (b), along the duct length. The averaged and fluctuating components of the streamwise velocity were normalized by the value of the average mass velocity of a single-phase flow at the duct inlet U m1 . Comparisons with the data of [39] were made for the 17th and 18th obstacles. The height of the duct with a square cross section was H = 60 mm. The profiles of the mean longitudinal velocity component agreed well with the experimental data (the difference did not exceed 5-7%). The agreement between the measurements and numerical predictions for longitudinal velocity pulsations was also quite good (the difference did not exceed 10%) except for the near-wall region. agreement between the measurements and numerical predictions for longitudinal velocity pulsations was also quite good (the difference did not exceed 10%) except for the nearwall region.
The results of measurements [40] and RANS numerical simulations with various isotropic turbulence models (k-ε, v2f, and k-ω shear stress tensor (SST)) [41] for the flow in the ribbed duct were used for heat transfer comparisons. Satisfactory agreement with the data of other authors for a single-phase flow around a two-dimensional obstacle was obtained (the maximum differences did not exceed 15%), except for the duct cross section near the upstream obstacle at x/h < 2 (see Figure 5). Here, Nu is the Nusselt number in a ribbed duct and Nu0 is the Nusselt number in a smooth duct for a single-phase flow. The Nusselt number at a constant value of heat flux density (qW = const) is determined by the formula: Comparisons with the data [40,41] were made for the 7th and 8th obstacles. All predictions were carried out for a flat duct with a square cross section and a height of H = 30 mm.    The results of measurements [40] and RANS numerical simulations with various isotropic turbulence models (k-ε, v2f, and k-ω shear stress tensor (SST)) [41] for the flow in the ribbed duct were used for heat transfer comparisons. Satisfactory agreement with the data of other authors for a single-phase flow around a two-dimensional obstacle was obtained (the maximum differences did not exceed 15%), except for the duct cross section near the upstream obstacle at x/h < 2 (see Figure 5). Here, Nu is the Nusselt number in a ribbed duct and Nu 0 is the Nusselt number in a smooth duct for a single-phase flow. The Nusselt number at a constant value of heat flux density (q W = const) is determined by the formula: (a) (b)   [40]; the lines are predictions: v2f, k and k-ε are predictions of [41], and RSM is the authors' simulations. Comparisons with the data [40,41] were made for the 7th and 8th obstacles. All predictions were carried out for a flat duct with a square cross section and a height of H = 30 mm.

The RANS Results and Discussion
All 2D numerical simulations were carried out for a mixture of air with water drops at the duct inlet for the case of a downward two-phase flow at atmospheric pressure. Ribs were installed on the "bottom" wall of the flat duct. All simulations were performed for the flow around the system of the 2nd and 3rd obstacles. The computational domain included two square ribs with a height of h = 4 mm. The height of a smooth duct was H = 40 mm (H/h = 10), and the distance between two ribs was p/h = 5-12. The massaverage gas velocity in the inlet cross section in the computational domain varied within U m1 = 5-20 m/s, and the Reynolds number for the gas phase, constructed from the massaverage gas velocity at the inlet and the duct height, was Re H = HU m1 /ν ≈ (0.6-5) × 10 4 . The initial average droplet diameter was d 1 = 5-50 µm, and their mass concentration was M L1 = 0-10%. The initial temperature of the gaseous and dispersed phases was A turbulent flow is 3D in nature. Nevertheless, there are many cases when it is possible to use a 2D approach to describe a quasi-two-dimensional turbulent flow, for example, if the duct width, Z, is much greater than its height, H (Z/H > 10). The authors of [42] recommended the consideration of the turbulent solid particle-laden flow in a backward-facing step in a flat channel as two-dimensional due to the large aspect ratio of Z/H.

Flow Structure
The streamlines for a gas-droplet flow around the system of two ribs are shown in Figure 6. The complex vortex structures of the averaged flow between two ribs are clearly visible. The formation of two regions of the flow recirculation is shown. The first large recirculation region formed behind the upstream rib due to the separation of the two-phase flow at the backward-facing step (BFS). A small corner vortex was located at the end of the reverse step. The second one formed due to the droplet-laden flow separation before the downstream rib (forward-facing step (FFS)) when the fluid flow left the cell between the two ribs. It was much shorter than the previous one. Figure 6. The complex vortex structures of the averaged flow between two ribs are clearly visible. The formation of two regions of the flow recirculation is shown. The first large recirculation region formed behind the upstream rib due to the separation of the twophase flow at the backward-facing step (BFS). A small corner vortex was located at the end of the reverse step. The second one formed due to the droplet-laden flow separation before the downstream rib (forward-facing step (FFS)) when the fluid flow left the cell between the two ribs. It was much shorter than the previous one.  The structure of a turbulent two-phase flow showed significant changes when flowing around a system of obstacles installed on one of the duct walls. The profiles of the averaged streamwise velocity components of the gaseous and dispersed phases were similar to those for the single-phase flow regime (see Figure 7a). The gas velocity in the gasdroplet flow was slightly (≤3%) ahead of the single-phase flow velocity. The drop velocity had the greatest value for the downward flow due to their inertia. Two regions with negative values for the longitudinal velocity of the gas-droplet flow are shown, which were confirmed by the data in Figure 6. The length of the main recirculation zone of the flow was xR1 ≈ 4.1h, and the length of the second recirculation region in front of the step ahead   Figure 7b shows the transverse distributions of the kinetic energy (TKE) of carrier phase turbulence for a 2D flow. The TKE was calculated by the formula for a two-dimensional flow: The highest turbulence values were obtained for the mixing layer. The level of kinetic energy of turbulence increased as the downstream obstacle approached. The maximum value of gas-phase turbulence was obtained at x/h = 9 (the upper corner of the downstream rib), and it was almost twice as high as the values for the TKE between the ribs. The turbulence of the flow was associated with the flow around the obstacle.
The dimensionless temperature distributions of the single-phase flow and the gas and dispersed phases are shown in Figure 7c. All profiles in Figure 7c are qualitatively similar to each other. The gas temperature in the gas-droplet flow was lower than the corresponding value for a single-phase flow due to droplet evaporation. Let us note that the droplet temperature profile for the first two sections, x/h = and 3, did not start from the wall (y/h = 0) as for the gas phase, but it is shifted from the wall by a small distance towards the flow core. This is explained by the absence of droplets in the near-wall zone in the area of flow separation due to their evaporation close to the wall between the ribs. The non-dimensional vorticity profiles are given in Figure 7d. They were calculated using the well-known formula: The magnitudes of vorticity were mainly negative values (because , except in the near-wall region inside the flow recirculation zone (see Figure 7d). The minimal values are shown in the outer shear layer of the separation zone and on the top wall of the downstream rib. The maximal positive value was observed close to the wall of the ribbed wall. In the case of two-phase mist flow, the magnitude of vorticity was slightly higher than that of the single-phase flow (up to 4%). Figure 8 shows the profiles of the dispersed-phase mass concentration, ML/ML1, for various droplet mass fractions (a) and their initial diameters (b). Obviously, due to the evaporation of droplets, their mass fraction decreased continuously, both streamwise and in traverse directions, when approaching the wall of the heated duct between the ribs. This was typical of the numerical data given in Figure 8a  The structure of a turbulent two-phase flow showed significant changes when flowing around a system of obstacles installed on one of the duct walls. The profiles of the averaged streamwise velocity components of the gaseous and dispersed phases were similar to those for the single-phase flow regime (see Figure 7a). The gas velocity in the gas-droplet flow was slightly (≤3%) ahead of the single-phase flow velocity. The drop velocity had the greatest value for the downward flow due to their inertia. Two regions with negative values for the longitudinal velocity of the gas-droplet flow are shown, which were confirmed by the data in Figure 6. The length of the main recirculation zone of the flow was x R1 ≈ 4.1h, and the length of the second recirculation region in front of the step ahead was x R2 ≈ 1.1h. The lengths of the recirculation zones were determined from the zero value of the flow velocity. Figure 7b shows the transverse distributions of the kinetic energy (TKE) of carrier phase turbulence for a 2D flow. The TKE was calculated by the formula for a twodimensional flow: The highest turbulence values were obtained for the mixing layer. The level of kinetic energy of turbulence increased as the downstream obstacle approached. The maximum value of gas-phase turbulence was obtained at x/h = 9 (the upper corner of the downstream rib), and it was almost twice as high as the values for the TKE between the ribs. The turbulence of the flow was associated with the flow around the obstacle.
The dimensionless temperature distributions of the single-phase flow and the gas and dispersed phases are shown in Figure 7c. All profiles in Figure 7c are qualitatively similar to each other. The gas temperature in the gas-droplet flow was lower than the corresponding value for a single-phase flow due to droplet evaporation. Let us note that the droplet temperature profile for the first two sections, x/h = and 3, did not start from the wall (y/h = 0) as for the gas phase, but it is shifted from the wall by a small distance towards the flow core. This is explained by the absence of droplets in the near-wall zone in the area of flow separation due to their evaporation close to the wall between the ribs. The non-dimensional vorticity profiles are given in Figure 7d. They were calculated using the well-known formula: The magnitudes of vorticity were mainly negative values (because ∂V ∂x ∂U ∂y ), except in the near-wall region inside the flow recirculation zone (see Figure 7d). The minimal values are shown in the outer shear layer of the separation zone and on the top wall of the downstream rib. The maximal positive value was observed close to the wall of the ribbed wall. In the case of two-phase mist flow, the magnitude of vorticity was slightly higher than that of the single-phase flow (up to 4%). Figure 8 shows the profiles of the dispersed-phase mass concentration, M L /M L1 , for various droplet mass fractions (a) and their initial diameters (b). Obviously, due to the evaporation of droplets, their mass fraction decreased continuously, both streamwise and in traverse directions, when approaching the wall of the heated duct between the ribs. This was typical of the numerical data given in Figure 8a,b. The distributions of the mass fraction of droplets with changes in their initial amounts had qualitatively similar forms (see Figure 8a). fraction of droplets with changes in their initial amounts had qualitatively similar forms (see Figure 8a).
A change in the initial diameter of the liquid droplets had a more complex effect on the course of the evaporation processes (see Figure 8b). In the flow core, this value trended toward the corresponding value at the inlet to the computational domain, and ML/ML1 → 1. This is explained by the almost complete absence of droplet evaporation. Fine particles at Stk < 1 penetrated into the region of flow separation and were observed over the entire cross section of the duct. Large inertial droplets (d1 = 100 μm, Stk > 1) almost did not penetrate into the flow recirculation zone, and they were present in the mixing layer and the flow core. In the near-wall zone, large drops were observed only behind the reattachment point. The largest and inertial droplets (d1 = 100 μm) accumulated in the near-wall region towards the downstream obstacle. Finely dispersed low-inertia droplets could leave the region between the two ribs due to their low inertia, while large drops could not leave this region. This led to an increase in the droplet mass fraction in this flow region and towards the downstream obstacle. In order to clearly display the flow structure in the inter-rib cavity, the contours of the nondimensional mean streamwise velocity, U/Um1 (a), and the temperature,  A change in the initial diameter of the liquid droplets had a more complex effect on the course of the evaporation processes (see Figure 8b). In the flow core, this value trended toward the corresponding value at the inlet to the computational domain, and M L /M L1 → 1. This is explained by the almost complete absence of droplet evaporation. Fine particles at Stk < 1 penetrated into the region of flow separation and were observed over the entire cross section of the duct. Large inertial droplets (d 1 = 100 µm, Stk > 1) almost did not penetrate into the flow recirculation zone, and they were present in the mixing layer and the flow core. In the near-wall zone, large drops were observed only behind the reattachment point. The largest and inertial droplets (d 1 = 100 µm) accumulated in the near-wall region towards the downstream obstacle. Finely dispersed low-inertia droplets could leave the region between the two ribs due to their low inertia, while large drops could not leave this region. This led to an increase in the droplet mass fraction in this flow region and towards the downstream obstacle.
In order to clearly display the flow structure in the inter-rib cavity, the contours of the nondimensional mean streamwise velocity, U/U m1 (a), and the temperature, , in two-phase mist flow are shown in Figure 9. Largescale and small-scale flow recirculation zones behind the upwind rib (BFS) and before the downstream rib (FFS) can be found in Figure 9a. The small corner vortex directly behind the upstream rib was observed. The length of the main recirculation zone of the flow was x R1 ≈ 4.1h, and the length of the second recirculation region in front of the step ahead was x R2 ≈ 1.1h. The lengths of the recirculation zones were determined from the zero value of the mean streamwise flow velocity (U = 0). In this region, the gas temperature increased, and it led to the suppression of heat transfer (see Figure 9b). These conclusions agree with the data of Figures 6 and 7a,c. Water 2022, 14, x FOR PEER REVIEW 13 of 19 ≈ 4.1h, and the length of the second recirculation region in front of the step ahead was xR2 ≈ 1.1h. The lengths of the recirculation zones were determined from the zero value of the mean streamwise flow velocity (U = 0). In this region, the gas temperature increased, and it led to the suppression of heat transfer (see Figure 9b). These conclusions agree with the data of Figures 6 and 7a,c.

Heat Transfer
The influence of the initial mass fraction (a) and droplet diameter (b) of the dispersed phase on the Nusselt number distribution in a two-phase flow along the duct length is shown in Figure 10. A significant HTE in the two-phase mist flow (up to 2.5 times) compared to a single-phase flow in a ribbed channel was obtained with the addition of evaporating water drops into a single-phase gas flow (see Figure 10a). Droplets of the minimum diameter (d1 = 5 µ m) evaporated most intensely, and the largest ones evaporated least intensely (d1 = 100 µ m) (see Figure 10b). The sizes of the zone of two-phase flow and the zone of HTE also decreased. This was an obvious fact for the evaporation of droplets in the two-phase mist flows, which was associated with a significant interface reduction;

Heat Transfer
The influence of the initial mass fraction (a) and droplet diameter (b) of the dispersed phase on the Nusselt number distribution in a two-phase flow along the duct length is shown in Figure 10. A significant HTE in the two-phase mist flow (up to 2.5 times) compared to a single-phase flow in a ribbed channel was obtained with the addition of evaporating water drops into a single-phase gas flow (see Figure 10a). Droplets of the minimum diameter (d 1 = 5 µm) evaporated most intensely, and the largest ones evaporated least intensely (d 1 = 100 µm) (see Figure 10b). The sizes of the zone of two-phase flow and the zone of HTE also decreased. This was an obvious fact for the evaporation of droplets in the two-phase mist flows, which was associated with a significant interface reduction; it was first shown by the authors of this work for a gas-droplet flow in a system of twodimensional obstacles. Heat transfer was attenuated and trended toward the corresponding value for the single-phase flow in the region of flow separation for the most inertial droplets. These drops did not penetrate into the flow separation region behind the upstream rib (BFS). An increase in heat transfer was obtained in the region behind the point of flow reattachment. A decrease in heat transfer was shown in the section of flow separation towards the downstream rib (FFS). The most inertial droplets also did not leave the region between the two ribs and accumulated in front of the downstream obstacle. it was first shown by the authors of this work for a gas-droplet flow in a system of twodimensional obstacles. Heat transfer was attenuated and trended toward the corresponding value for the single-phase flow in the region of flow separation for the most inertial droplets. These drops did not penetrate into the flow separation region behind the upstream rib (BFS). An increase in heat transfer was obtained in the region behind the point of flow reattachment. A decrease in heat transfer was shown in the section of flow separation towards the downstream rib (FFS). The most inertial droplets also did not leave the region between the two ribs and accumulated in front of the downstream obstacle. The effect of the gas Reynolds number, Re, and the initial mass fraction of the dispersed phase, ML1, on the thermal hydraulic performance parameter is shown in Figure  11. The wall friction coefficient, Cf, was calculated using the formula ( ) Here, Nu0 and Cf0 are the maximal Nusselt number and wall friction coefficient in the two-phase mist flow of a fully developed smooth duct, other conditions being equal. Nu/Nu0/(Cf/Cf0) is the thermal hydraulic performance parameter. This is the ratio of the maximal Nusselt numbers divided by the maximal wall friction coefficient ratio. The ribbed surface provided a much better thermohydraulic performance than a smooth duct in the case of a droplet-laden turbulent mist flow, with other conditions being identical. This effect was quite pronounced at small Reynolds number values of Re < 10 4 . It should be noted that the wall friction coefficient ratio, Cf/Cf0, was taken to the first power.  The effect of the gas Reynolds number, Re, and the initial mass fraction of the dispersed phase, M L1 , on the thermal hydraulic performance parameter is shown in Figure 11. The wall friction coefficient, C f , was calculated using the formula C f /2 = τ W / ρU 2 m1 . Here, Nu 0 and C f0 are the maximal Nusselt number and wall friction coefficient in the two-phase mist flow of a fully developed smooth duct, other conditions being equal. Nu/Nu 0 /(C f /C f0 ) is the thermal hydraulic performance parameter. This is the ratio of the maximal Nusselt numbers divided by the maximal wall friction coefficient ratio. The ribbed surface provided a much better thermohydraulic performance than a smooth duct in the case of a dropletladen turbulent mist flow, with other conditions being identical. This effect was quite pronounced at small Reynolds number values of Re < 10 4 . It should be noted that the wall friction coefficient ratio, C f /C f0 , was taken to the first power. The effect of the gas Reynolds number, Re, and the initial mass fraction of the dispersed phase, ML1, on the thermal hydraulic performance parameter is shown in Figure  11. The wall friction coefficient, Cf, was calculated using the formula ( ) Here, Nu0 and Cf0 are the maximal Nusselt number and wall friction coefficient in the two-phase mist flow of a fully developed smooth duct, other conditions being equal. Nu/Nu0/(Cf/Cf0) is the thermal hydraulic performance parameter. This is the ratio of the maximal Nusselt numbers divided by the maximal wall friction coefficient ratio. The ribbed surface provided a much better thermohydraulic performance than a smooth duct in the case of a droplet-laden turbulent mist flow, with other conditions being identical. This effect was quite pronounced at small Reynolds number values of Re < 10 4 . It should be noted that the wall friction coefficient ratio, Cf/Cf0, was taken to the first power.

Comparison with Results of Other Authors
Comparisons with the data of LES simulations of a solid particle-laden flow around a two-dimensional obstacle were made according to the conditions of [10]. The following data were used for the comparative analysis: h = H/7, V P1 = U m1 /25, the Reynolds number was plotted from the obstacle height, Re = HU m1 /ν = 740, ρ P /ρ = 769.2, and ρ P was the particle material density. The height of the boundary layer for a single-phase flow in the inlet section of the computational domain was δ = 7h, and the carrier phase was atmospheric air at T = 293 K (see Figure 12). Here, h = 7 mm was the obstacle height, H = 1 mm, U m1 = 1.59 m/s was the free flow velocity, and V P1 = 0.06 m/s. The two-dimensional obstacle was square in cross section and was mounted on the bottom wall. The flow of solid particles was blown vertically through a flat slot along the normal surface at distance h from the trailing edge of the obstacle. The number of solid particles during the LES calculation was 2 × 10 5 . The calculations were performed for three Stokes numbers, St + = τu * 2 /ν = 0.25, 1, 5, and 25, where τ = ρ P d 2 /(18µ) was the particle relaxation time and u * = 0.5 m/s was the friction velocity for a single-phase flow without particles, other things being equal. This corresponded to the solid particle diameters d = 8, 15, 34, and 76 µm. The calculations were carried out in a two-dimensional formulation for an isothermal two-phase flow around a single obstacle.
ticle material density. The height of the boundary la section of the computational domain was δ = 7h, a air at T = 293 K (see Figure 12). Here, h = 7 mm wa 1.59 m/s was the free flow velocity, and VP1 = 0.06 m square in cross section and was mounted on the bo was blown vertically through a flat slot along the trailing edge of the obstacle. The number of solid pa 2 × 10 5 . The calculations were performed for three S and 25, where τ = ρPd 2 /(18µ) was the particle relaxa tion velocity for a single-phase flow without particl responded to the solid particle diameters d = 8, 15, carried out in a two-dimensional formulation for an single obstacle. The profiles of the dispersed phase concentrat at y = 0.02h are shown in Figure 13. Here, Cb is the the hole (slot) width at the inlet to the computation increased, heavier particles stopped penetrating int lower concentrations along the obstacle wall. The lo for all studied Stokes numbers (particle diameters particle flow. A characteristic feature of the low-ine in the concentration of particles near the obstacle w [10][11][12][13][14][15][16][17][18][19][20]. Most likely, such an accumulation of partic obstacle can be explained by the effect of the accu Eulerian simulations, an increase in concentration w The profiles of the dispersed phase concentration in the near-wall zone of the plate at y = 0.02h are shown in Figure 13. Here, C b is the mean concentration of particles over the hole (slot) width at the inlet to the computational domain. As the Stokes number, St + , increased, heavier particles stopped penetrating into the recirculation region, resulting in lower concentrations along the obstacle wall. The local maximum concentration at x/h ≈ 1 for all studied Stokes numbers (particle diameters) is explained by the injection of the particle flow. A characteristic feature of the low-inertia particles was a significant increase in the concentration of particles near the obstacle wall, according to the LES data (C/C b ≈ [10][11][12][13][14][15][16][17][18][19][20]. Most likely, such an accumulation of particles in the corner near the wall of the obstacle can be explained by the effect of the accumulation of particles in [42]. For our Eulerian simulations, an increase in concentration was also obtained, but the values were much smaller (by a factor of approximately 8-10). For inertial particles at St + = 5, the region turned out to be almost completely free of solid particles. This was typical for both the data of the LES calculations [10] and our numerical calculations. Behind the obstacle, a decrease in the particle concentration in the near-wall region was observed, and here our numerical predictions agreed satisfactorily with the LES data (the difference did not exceed 20% at St + = 1 and 5 and did not exceed 100% at St + = 0.25).  Figure 13 shows the concentration profiles of the dispersed phase when the Stokes number, St+, is varied along the length of the channel behind a two-dimensional obstacle. Particles at St+ = 0.25 accumulated in the near-wall region near the bottom wall. Further downstream, heavier particles gradually left the recirculation region, and at St+ > 1 the decrease in their distribution profile was similar to a Gaussian distribution. For the largest particles at St+ = 25, according to the results of our numerical predictions, an underestimation of the position of the concentration maximum was observed, and in general the particles rose lower than according to the LES results [10]. Figure 14 shows the profiles of the dispersed phase concentration when the Stokes number, St+, was varied along the length of the duct behind a two-dimensional obstacle. Particles at St+ = 0.25 accumulated in the near-wall region near the bottom wall. Further downstream, heavier particles gradually left the recirculation region, and at St+ > 1 the decrease in their distribution profile was similar to a Gaussian distribution. An underestimation of the position of the concentration maximum was observed, according to the results of our numerical predictions for the largest particles at St+ = 25. The maximal penetration coordinate in the transverse directions in our RANS predictions was smaller than that in the LES results [10]. Figure 14. The transverse profiles of particle concentrations for various Stokes numbers, St+, after a 2D obstacle. The points are LES calculations [10]; the lines are the authors' predictions.

Conclusions
Two-dimensional numerical simulations of the local flow structure, turbulence, and heat transfer in a ribbed flat duct during the evaporation of water droplets in a gas flow  Figure 13 shows the concentration profiles of the dispersed phase when the Stokes number, St + , is varied along the length of the channel behind a two-dimensional obstacle. Particles at St + = 0.25 accumulated in the near-wall region near the bottom wall. Further downstream, heavier particles gradually left the recirculation region, and at St + > 1 the decrease in their distribution profile was similar to a Gaussian distribution. For the largest particles at St + = 25, according to the results of our numerical predictions, an underestimation of the position of the concentration maximum was observed, and in general the particles rose lower than according to the LES results [10]. Figure 14 shows the profiles of the dispersed phase concentration when the Stokes number, St + , was varied along the length of the duct behind a two-dimensional obstacle. Particles at St + = 0.25 accumulated in the near-wall region near the bottom wall. Further downstream, heavier particles gradually left the recirculation region, and at St + > 1 the decrease in their distribution profile was similar to a Gaussian distribution. An underestimation of the position of the concentration maximum was observed, according to the results of our numerical predictions for the largest particles at St + = 25. The maximal penetration coordinate in the transverse directions in our RANS predictions was smaller than that in the LES results [10].  Figure 13 shows the concentration profiles of the dispersed phase when the Stokes number, St+, is varied along the length of the channel behind a two-dimensional obstacle. Particles at St+ = 0.25 accumulated in the near-wall region near the bottom wall. Further downstream, heavier particles gradually left the recirculation region, and at St+ > 1 the decrease in their distribution profile was similar to a Gaussian distribution. For the largest particles at St+ = 25, according to the results of our numerical predictions, an underestimation of the position of the concentration maximum was observed, and in general the particles rose lower than according to the LES results [10]. Figure 14 shows the profiles of the dispersed phase concentration when the Stokes number, St+, was varied along the length of the duct behind a two-dimensional obstacle. Particles at St+ = 0.25 accumulated in the near-wall region near the bottom wall. Further downstream, heavier particles gradually left the recirculation region, and at St+ > 1 the decrease in their distribution profile was similar to a Gaussian distribution. An underestimation of the position of the concentration maximum was observed, according to the results of our numerical predictions for the largest particles at St+ = 25. The maximal penetration coordinate in the transverse directions in our RANS predictions was smaller than that in the LES results [10]. Figure 14. The transverse profiles of particle concentrations for various Stokes numbers, St+, after a 2D obstacle. The points are LES calculations [10]; the lines are the authors' predictions.

Conclusions
Two-dimensional numerical simulations of the local flow structure, turbulence, and heat transfer in a ribbed flat duct during the evaporation of water droplets in a gas flow Figure 14. The transverse profiles of particle concentrations for various Stokes numbers, St + , after a 2D obstacle. The points are LES calculations [10]; the lines are the authors' predictions.

Conclusions
Two-dimensional numerical simulations of the local flow structure, turbulence, and heat transfer in a ribbed flat duct during the evaporation of water droplets in a gas flow were carried out. The set of steady-state RANS equations written with consideration of the influence of droplet evaporation on the transport processes in gas is used. The two-fluid Eulerian approach was used to describe the flow dynamics and heat and mass transfer in the dispersed phase. To describe the turbulence of the gas phase, an elliptical blending RSM model was employed.
It was shown that the transverse profiles of the averaged longitudinal velocity components of the gaseous and dispersed phases were similar to those of the single-phase flow regime. The gas velocity in the gas-droplet flow was slightly (≤3%) higher than those in the single-phase flow. The droplet velocity is higher than the gas-phase velocity in the two-phase flow. Finely dispersed droplets (Stk < 1) penetrated well into the region of flow recirculation and were observed over the entire cross section of the duct. They could leave the region between the two ribs due to their low inertia. Large inertial droplets (Stk > 1) were present only in the mixing layer and the flow core and accumulated in the near-wall region close to the downstream wall of the rib. A significant increase in heat transfer (up to 2.5 times) in comparison with a single-phase flow in a ribbed duct was shown when evaporating water drops were added to a single-phase gas turbulent flow. For the most inertial droplets, which were not involved in the separation motion in the region of the main recirculation zone behind the BFS (upstream rib), the heat transfer intensification decreased and trended toward the corresponding value for the single-phase flow regime in the recirculation zone. An increase in heat transfer was obtained behind the reattachment point. A decrease in heat transfer was shown in the zone close to the FFS (downstream rib).