**Experimental and Numerical Modeling of Fluid Flow**

Editor

**R ¨udiger Schwarze**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editor* Rudiger Schwarze ¨ TU Bergakademie Freiberg Germany

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Applied Sciences* (ISSN 2076-3417) (available at: https://www.mdpi.com/journal/applsci/special issues/Experimental Numerical Modeling Fluid Flow).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-5469-3 (Hbk) ISBN 978-3-0365-5470-9 (PDF)**

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

### **Contents**


### *Editorial* **Experimental and Numerical Modeling of Fluid Flow**

**Rüdiger Schwarze**

Institute of Mechanics and Fluid Dynamics, Lampadiusstraße 4, 09596 Freiberg, Germany; ruediger.schwarze@imfd.tu-freiberg.de; Tel.: +49-3731-392486

### **1. Introduction**

Fluid dynamics is often related to complex flow conditions and systems, either in the context of fundamental research or in the context of industrial processes. Typically, the study of such flows requires experimental or numerical models which capture their complex, i.e., multiphase and turbulent, nature. Although such models have been the state of the art for a long time in areas such as aerodynamics or weather simulation, the field of the application of experimental and numerical flow models is constantly expanding. The goal of this Special Issue is to present an overview of applied experimental and numerical flow models, with a strong focus on the following topics: (i) new areas of application, (ii) advanced experimental or numerical models, (iii) innovative modeling approaches, and (iv) challenges in modeling techniques.

### **2. New Applications of Experimental and Numerical Flow Models**

T. Zisis, K. Vasilopoulos, and I. Sarris [1] investigate the prevailing conditions in a railway tunnel after a train fire accident by means of CFD simulations. The goal is to analyze the ability of the ventilation system to create the proper conditions for safe passenger evacuation. Three scenarios with different fan activation times and different evacuation processes were examined. They found that the most important action in a tunnel fire is the time at which the ventilation system is activated after the start of a fire.

C. L. Pavlidis, A. V. Palampigik, K. Vasilopoulos, I. C. Lekakis, and I. E. Sarris [2] examine the airflow and pollutant dispersion around an isolated cubical building located in a warm Mediterranean climate with CFD simulations adapted to the local micro-climate conditions. The performance of the numerical model is tested with the help of measured data from the SILSOE cube experiment. In the simulations, the influence of the thermal boundary conditions of the building on the airflow and the pollutant dispersion is evaluated.

M. Zaboli, S. S. Mousavi Ajarostaghi, and M. Saffari Pour [3] analyze different configurations of a parabolic trough solar collector with inner helical axial fins as a swirl generator. After a successful validation of their CFD model, they use CFD simulations in order to resolve the inner flow and temperature fields. They show that the thermal performance of the collector can be significantly improved by using an innovative collector design.

K. Yousef, A. Hegazy, and A. Engeda [4] study complex flow behavior during gas entrainment in water in an inverted vertical U-tube with the addition of side air or water vapor. The results from both experiments and CFD simulations reveal a complex interplay between the water mass flow rate, pressure, and extracted air or vapor flow rates. Furthermore, critical states of the system are identified.

Z. Zhu, Y. Ju, and C. Zhang [5] investigate the obstructive sleep apnea–hypopnea syndrome, a highly prevalent respiratory disorder, in an experimental approach based on a computed tomography scanned extra-thoracic airway model. They use particle image velocimetry to study the air flow in their model. The in vitro measurements of oscillatory respiratory flow velocity show the temporal evolution of the complex flow patterns and

**Citation:** Schwarze, R. Experimental and Numerical Modeling of Fluid Flow. *Appl. Sci.* **2022**, *12*, 9042. https://doi.org/10.3390/app12189042

Received: 5 September 2022 Accepted: 7 September 2022 Published: 8 September 2022

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

**Copyright:** © 2022 by the author. 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 potential wall collapse of the extra-thoracic airway model with obstructive sleep apnea– hypopnea syndrome.

D. K. Kim [6] analyzes the performance improvement of refrigerators by different innovative methods. The effects of three methods, modifying the flow channel and the operating conditions and changing the positions of the fan in the mechanical chamber room, are analyzed in a combined experimental and numerical approach. The results show that the efficiency of the refrigerator can be improved by these methods.

T. Bano, F. Hegner, M. Heinrich, and R. Schwarze [7] study the fluid–structure interaction (FSI) involving flexible elastic structures in a water cross stream with an FSI-CFD model. Transient simulations are carried out for flexible flaps of different thicknesses in a laminar glycerin flow with Reynolds numbers ranging from 3 < Re < 12. The bending line and the maximum tip deflection of the flaps in the simulations agree well with findings from a previous experimental study.

F. Wang, L. Wang, G. Chen, and D. Zhu [8] examine the oil droplet size distribution in an aero-engine bearing chamber with CFD and population balance models. The CFD simulations reveal the correlations between the initial oil droplet size distribution, the oil droplets' coalescence and breakup, air, and oil mass flow, analyzed using the CFD and population balance models.

X. Qu, Q. Guo, Y. Zhang, X. Qi, and L. Liu [9] propose a miniaturized four-sensor electrical probe with a new signal processing method as an innovative multiphase flow measurement technique. Measurements were taken during a gas–liquid two-phase flow experiment in cap-bubble flow regime to test the performance of the proposed signal process scheme. The local flow parameters obtained by the probe measurement are in good agreement with the results from the visual measurement technique in the same flow conditions, proving the correctness of the proposed method.

Finally, Y.-B. Seo, S. S. Paramanantham, J.-Y. Bak, B. Yun, and W.-G. Park [10] use a CFD approach to investigate a boiling model for rod bundle flows related to reactor cooling in nuclear power plants. The comparison of the numerical results with the experimental values of the vapor volume fraction and vapor bubble–water interfacial area concentration show good agreement for different initial conditions.

**Funding:** This research received no external funding.

**Acknowledgments:** Thanks to all authors and peer reviewers for their valuable contributions to this Special Issue. I would also like to express my gratitude to all the staff and people involved in this Special Issue.

**Conflicts of Interest:** The author declares no conflict of interest.

### **References**


### *Article* **Numerical Simulation of a Fire Accident in a Longitudinally Ventilated Railway Tunnel and Tenability Analysis**

**Thomas Zisis, Konstantinos Vasilopoulos \* and Ioannis Sarris**

Laboratory of Thermo Fluid Systems (LTFS), Department of Mechanical Engineering, Ancient Olive Grove Campus, University of West Attica, Thivon Str. 250, Egaleo, 12244 Athens, Greece; msrtf20x03@uniwa.gr (T.Z.); sarris@uniwa.gr (I.S.)

**\*** Correspondence: kvassil@uniwa.gr; Tel.: +30-2105381131

**Abstract:** The present study examines the prevailing conditions in a railway tunnel after a train fire accident and the ability of the ventilation system to create the proper conditions for a safe passenger evacuation. The examined scenario included an event of a 20-MW diesel pool fire on a suburban train, immobilized in the middle of a 1.5-Km long, linear shaped rectangular tunnel ventilated by a longitudinal jet fan system, and the people's movement during the evacuation was effectuated along walking platforms. More specifically, three scenarios with different fan activation times and different evacuation processes were examined. A Large Eddy simulation model (LES) was used for the simulation of the air flow in the railway tunnel. The evaluation of the ventilation system criteria considered the achievement of the air critical velocity inside the railway tunnel, and for the people's safe evacuation, the Fractional Effective Dose (FED) value was examined. It was found that the most important action in a tunnel fire is the time, after the start of a fire, the ventilation system is activated.

**Keywords:** railway tunnel fire; FDS; tunnel ventilation; tenability analysis; numerical simulation

### **1. Introduction**

Tunnels are complex infrastructures, usually of semicircular or rectangular crosssections, constructed to establish connections between different sections of a road or railroad network. The growing reports of indoor fire accidents with high population mobility, such as tunnels, have created the need to consider the prevailing conditions in such cases as well as the best evacuation plan to ensure human life [1].

The consequences after a tunnel fire can be much more serious than an outdoor fire because of high smoke concentrations and confined movement space [2,3]. Accidents involving fire and hazardous material displacement can be studied for different urban scales and settings. Although the fire poses a threat near its vicinity, toxic substances can travel along the flow of the air [4]. After a fire accident, high temperatures and toxic fumes are produced, which have harmful consequences for passengers who are trapped inside the tunnel. In 2003, an arson fire in the subway tunnel in Daegu, South Korea resulted in 198 deaths. In 2017, the Jinji Road Tunnel Fire in Shanxi, China, 2017, caused 31 deaths [5]. A wrong evacuation plan can lead people in to the fire or smoke direction, with a negative effect on their health due to the presence of toxic gases, heat, and radiation [6,7].

To minimize human exposure to smoke during a tunnel evacuation, various types of ventilation systems are used to remove the smoke from the tunnel or to keep the smoke in the tunnel above the head [8,9]. These can be either mechanical or natural ventilation systems [10]. Forced ventilation systems use fans to extract or push smoke out of the tunnel. Natural ventilation with vertical wells in the roof of the tunnel take advantage of the buoyancy of hot smoke, which escapes through the shafts due to the stack effect and piston effect due to the movement of vehicles or trains inside the tunnel. This approach generally does not require mechanical fans, resulting in space-saving at the higher level of the tunnel, and it is suitable for overground [11].

**Citation:** Zisis, T.; Vasilopoulos, K.; Sarris, I. Numerical Simulation of a Fire Accident in a Longitudinally Ventilated Railway Tunnel and Tenability Analysis. *Appl. Sci.* **2022**, *12*, 5667. https://doi.org/10.3390/ app12115667

Academic Editor: Rüdiger Schwarze

Received: 5 May 2022 Accepted: 31 May 2022 Published: 2 June 2022

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

**Copyright:** © 2022 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/).

One of the key criteria for the passenger's safety after a fire accident is the number of toxic products of combustion that the individual inhales during the event and the evacuation process. This is calculated using the FED (Fractional Effective Doses) index. The FED index was first introduced by Purser [12] and has since been adopted by various organizations, such as the NFPA [13]. The basic idea is that a person who inhales a certain dose of toxic products for a period of time is considered unconscious or dead. To assess how a person's life is threatened, FED is the main factor considered [1], alongside the thermal exposure and the low visibility due to smoke (hence the difficulty of accessing emergency exits).

Evacuation in underground public transport systems can be complicated and the outcome of an evacuation scenario depends on many different parameters. These parameters can be the people's walking speed, the visibility, the design of the train and tunnel, or the physical condition of the people [14].

Tunnels are designed with elevated platforms to assist a safe passenger evacuation in case of a fire accident. The reduced height difference between the level of the train and the platform, compared to the evacuation of a train directly to the rail area, is a method of increasing the safety of the population. The difficulty of evacuating a train to the rails, especially for people with mobility difficulties, has been analyzed in the METRO project of Ingason et al. [15]

The purpose of a tunnel ventilation system, besides for maintaining air quality, is to keep the longest part of a tunnel free of smoke for as long as possible. In the case of a fire accident inside a tunnel, passengers are exposed to smoke, even though the ventilation system exists. In such a situation, the activation of the ventilation system must be delayed slightly or operated at a slower speed so as not to disturb the natural stratification of the air and smoke layers [16,17]. This state will be maintained if the air velocity does not exceed the critical velocity [18]. The critical velocity is defined as the minimum longitudinal ventilation velocity to avoid the phenomenon of smoke reversal or backlayering in case of fire [19]. Depending on the air velocity, hot smoke and cold air are mixed, and thus the smoke is dispersed downstream of the fire throughout the tunnel [20]. Hot smoke also moves against the flow of air, upstream of the fire front, and as it cools it descends to the ground. This phenomenon is called smoke backlayering. The amount of smoke that will be inverted upstream of the fire front and the downstream distance of where it will remain depends on the ventilation conditions [19,21].

Hu et al. [22] suggested that smoke stops spreading when the horizontal inertia forces caused by the buoyancy of hot smoke are equal to the inertia forces of the ventilation air. Based on a series of small-scale experiments that studied the correlation between critical velocity and heat release rate (HRR—Heat Release Rate), for different geometries of fire sources, conducted by Oka and Atkinson [23] and Li et al. [24], they found that the dimensionless critical velocity is proportional to the dimensionless HRR1/3 in the case of low and medium HRR and is independent in the case of large HRR. Riess and Bettelini [25] found that the stack effect should be considered in the design of the ventilation of sloped tunnels with angles of inclination of more than 1–2%. Wu and Bakar [26] observed that, with the same fuel and longitudinal ventilation speed, the length of the inverted smoke layer gradually increases as the tunnel cross-section increases when the tunnel cross-section reaches a certain value.

Lee and Ryou [27] investigated the critical velocities for various tunnel aspect ratios in a series of small-scale experiments, and they found that the critical velocity increases with the aspect ratio. The study of a fire accident with vehicles blocking the tunnel showed [28,29] that, as the distance from the blockage upstream of the source fire accident increases, the length of the inverted smoke layer and the critical velocity initially decreases until it reaches a certain value that is equal to the value that these two sizes would have without the presence of an obstacle. Tanaka et al. [30] conducted a series of small-scale experiments to study the characteristics of a fire's, smoke considering the escalation of heat conduction through the tunnel walls and presented some equations for estimating the critical velocity and length of the inverted smoke layer.

According to NFPA [13], the critical velocity is "the minimum steady-state velocity of the ventilation airflow towards the fire, in a tunnel or corridor, required to prevent smoke reversal in the fire area." In this study, three different scenarios (A, B, and C), regarding the jet fan activation time, were examined. The jet fans were activated simultaneously at t = +180 s, +300 s, and +480 s, after the fire ignition, respectively. Then, the calculated volume fractions of toxic gases from the numerical simulation were used to calculate the FED value of the passengers evacuating the train. Three different scenarios for the passengers' evacuation were examined, with evacuation initiation times being t = +60 s, +120 s, and +300 s after the fire ignition.

Section 2 describes the numerical model, the dimensions and boundary conditions of the simulation model geometry, the grid formulation, and the grid independence test. Section 3 describes the validation method used to validate the numerical simulation results using empirical equations that resulted from experiments. The tenability analysis takes place in Section 4, where the mean and max FED values are compared for each scenario.

The literature, in its majority, studies the evacuation times, choke points during the evacuation of subway systems and the comfort of passengers evacuating tunnels via walking platforms, or how the critical velocity is affected by different conditions [14,24,29,31,32]. The present study focused on the correlation between the two actions taken in a tunnel fire accident: ventilation system activation and evacuation process initialization. These two actions play a great role in the exposure of the population to harmful gases and hot temperatures. The results show that the most important action is the activation of the ventilation system as soon as possible from the ignition of the fire. When this is not possible, according to the results, the best strategy is to wait for the people to evacuate the tunnel and then activate the ventilation system to prevent passengers from inhaling smoke due to the entertainment of the smoke layer in the tunnel ceiling.

### **2. Numerical Model**

The Fire Dynamics Simulation (FDS) is a code developed by the NIST in cooperation with the VIT research center of Finland. The code numerically solves the Navier-Stokes equations and finds application in low velocity thermally driven flows, smoke production, and mass and heat transfer. The FDS code is a low Mach number Large Eddy Simulation solver (the Mach number is under 0.3). To define the turbulent characteristics, the Smagorinsky model is applied.

The governing equations for mass and momentum conservation are:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \rho \stackrel{\rightarrow}{u} = \dot{m}\_b''' \tag{1}$$

$$\frac{\partial}{\partial t}(\rho \stackrel{\rightarrow}{u}) + \nabla \cdot \rho \stackrel{\rightarrow}{u}\stackrel{\rightarrow}{u} + \nabla \cdot \stackrel{\rightarrow}{p} = \rho \lnot + \stackrel{\rightarrow}{f}\_b + \nabla \cdot \pi\_{ij} \tag{2}$$

where *ρ* is the density, *g* is the gravity acceleration, *u* is the velocity, *p* is the pressure, *t* is the time, and *τij* = *μ* - <sup>2</sup>*Sij* <sup>−</sup> <sup>2</sup> <sup>3</sup> *δij*(∇*u*) is the stress tensor, where *μ* is the dynamic viscosity of the fluid, *Sij* is the strain tensor, and *δij* is the Kronecker delta.

The coefficient of turbulent consistency in Smagorinsky's analysis is defined as:

$$\mu\_{LES} = \rho(\mathcal{C}\_{s}\Lambda) \left( 2\mathcal{S}\_{i\bar{j}} \, : \, \mathcal{S}\_{i\bar{j}} - \frac{2}{3} (\nabla \overline{\mathbf{u}})^{2} \right)^{\frac{1}{2}} \tag{3}$$

where Δ = (*δxδyδz*) 1 <sup>3</sup> is the dimension of the spatial filter proportional to the dimension of the grid.

The energy conservation equation is written as:

$$\frac{\partial}{\partial t}(\rho h\_s) + \nabla \rho h\_s u = \frac{Dp}{Dt} + q''' - \dot{q}\_{b}''' - \nabla q'' + \varepsilon \tag{4}$$

where *hs* is the sensible enthalpy, . *q <sup>b</sup>* is the energy transferred to particles during evaporation (change of phase), and *q* are the radiative and conductive heat fluxes.

The equation for species conservation is written as:

$$\frac{\partial}{\partial t}(\rho \mathbf{Y}\_a) + \nabla \cdot \rho \mathbf{Y}\_a \overset{\rightarrow}{\boldsymbol{u}} = \nabla \rho \mathbf{D}\_a \nabla \mathbf{Y}\_a + \dot{m}\_a^{\prime\prime} + \dot{m}\_b^{\prime\prime} \tag{5}$$

where . *m <sup>b</sup>* = ∑ *a* . *m <sup>b</sup>*,*<sup>a</sup>* is the species production rate as particles or droplets.

The time step used is variable and adjusted by the FDS code. During the calculation, the time step is adjusted so that the CFL (Courant, Friedrichs, Lewy) condition is satisfied. The default value of the size of the time step is [31]:

$$DT = \frac{5(\delta \mathbf{x} \delta y \delta z)^{1/3}}{\sqrt{gH}} \tag{6}$$

where *δx*, *δy*, and *δz* are the dimensions of the smallest mesh cell, *H* is the height of the computational domain, and *g* is the acceleration of gravity.

### *2.1. Model Details*

The studied tunnel is made of rectangular cross-sections with dimensions typical for a double railway tunnel (y = 11 m, z = 7 m), and its length is 1.5 Km.

In all cases, the center of the fire source was located in the middle section of the tunnel, as shown in Figure 1 at (x = 795.5 m, y = 2.25 m, z = 4 m), and the intensity was 20 MW, which corresponds to a traction engine fire. The fire's surface was 1 m<sup>2</sup> and it was located on the train's roof. The surface of the fire was modeled as a burner type and was inserted into the model as an open vent. The fire growth followed a typical *at*<sup>2</sup> rule, where the growth rate *a* was 0.1876 kW/s<sup>2</sup> and achieved nominal HRR within 327 s. C12H23 was used as the fuel, which is the standard composition of diesel, with a yield of combustion products of 0.1 for CO and 0.09 for soot. Combustion efficiency is defined as the amount produced per gram of fuel diffused along with the flow, compared to the diffusion of fire into the environment.

**Figure 1.** Side view of the simulation tunnel, tunnel dimensions, jet fans, train, and fire source location.

The jet fans used in the simulations are 3 m long with a surface area of 1 m2. Four pairs of them were placed every 300 m at the height of 5 m and 1 m off the side walls

According to McGrattan et al. [32], the non-dimensional expression *D\**/*δx* was used to assess mesh resolution, where *δx* is the nominal size grid (in m), and the characteristic length of *D\** is expressed as follows

$$D^\* = \left(\frac{\dot{Q}}{\mathbb{C}\_p T\_{\infty} \rho\_{\infty} \sqrt{\mathcal{g}}}\right)^{2/5} \tag{7}$$

where *Q* is the fire heat release rate (HRR) (W), *ρ*<sup>∞</sup> is the ambient density (kg/m3), *Cp* is the specific heat capacity (J/kgK), *T*∞ is the ambient temperature (K), and *g* is the gravitational acceleration constant (m/s2). In the NUREG-1824 Directive published by the United States Nuclear Energy Regulatory Commission (USNRC) and the Electricity Research Institute (EPRI), it defines that *D*\* values should be in a range between 4–16 for an optimal analysis in FDS Code [33]. Unlike the number of cells, the grid resolution does not depend on the volume of the building, so the corresponding standard has the advantage of universal application regardless of the size of the building.

The grid independence test was conducted for three different meshes: coarse mesh with a mesh size of 0.75 m, medium mesh with a mesh size of 0.5 m, and fine mesh with a mesh resolution of 0.25 m. The number of cells for each type of mesh was 269,865, 924,000, and 4,254,149, respectively. As shown in Figure 2, the fine and medium meshes were compared. The grid errors using two computational grids were calculated using the following equation [34]

$$GCI = \frac{f\_2 - f\_1}{1 - r^p} \tag{8}$$

where *f*<sup>2</sup> is the numerical solution that results from the medium computational grid, and the corresponding from the finer is *f*1, *r* is the refinement factor between the two computational grids, and *p* is the accuracy of the algorithm, which is 3 for the present study.

**Figure 2.** (**a**) Comparison of the grid convergence index (GCI) for the U mean velocity and (**b**) mean carbon monoxide concentration for the medium and fine mesh in the symmetry plane of the tunnel, with coordinates (x = 1100 m, y = 5.5 m, and z = from 0 to 7 m) and 28 sampling points.

The grid independence examined the U velocity and carbon monoxide concentration values. The values were extracted by a line, downwind of the fire source, in the symmetry plane of the tunnel with coordinates (x = 1100 m, y = 5.5 m, and z = from 0 to 7 m) with 28 sampling points.

The results did not show great deviation, so a combination of the fine and medium mesh was used as the final mesh for the examination of the three main scenarios. The fine mesh was used for the field around the fire, while, for the rest of the tunnel, the medium mesh was used. For the sake of convenience, a mesh size of 0.25 m was used in the area around the fire (50 m downstream and 50 m upstream), while, for the rest of the tunnel, a mesh size of 0.5 m was used, as shown in Figure 3. The total number of cells was 1,355,200.

**Figure 3.** Grid meshes and fire source. The medium grid of 0.5 m was used for the majority of the tunnel, while the fine mesh of 0.25 was used near the field of the fire.

The mean value of errors of the U velocity for the coarse and medium meshes was approximately 14%, for the medium and fine mesh it was 5%, and for the medium and final mesh used it was 4%. For the CO concentrations, the mean error values for the coarse and medium mesh was 16%, for the medium and fine mesh was 8%, and for the medium and final mesh was 6% The mesh resolution, which was proposed in the UNSRC and EPRI directive, was used in many other tunnel fire simulations such as [35–37].

The work for the numerical simulation was conducted in the National HPC facility— ARIS, with computational time granted from the National Infrastructures for Research and Technology S.A. (GRNET S.A.). For each scenario, 22 cores and 56 GB of memory were used, and the typical time taken to complete each computation was 24 h.

### *2.2. Validation of Methodology*

Kurioka et al. [38] performed experiments using three different tunnel scales, 1:10, 1:2, and a full scale, with rectangular and horseshoe cross-sections. The fire sources used were of square cross section. The aspect ratio of the tunnel cross-section dimensions, HRR, and forced longitudinal ventilation velocity was variable. From the experimental data of the 1/10 scale model, empirical formulas were exported to calculate the flame inclination, the flame height, the maximum temperature of the smoke layer, and its position. The effect of tunnel cross-section dimensions was incorporated into these models. The ratio H3/2/b1/2, where H is the height of the tunnel model and b is the width, was found to be representative for the investigation of fire phenomena in tunnels. It was also confirmed that these empirical formulas were sufficiently applicable to predict phenomena occurring during a fire accident, close to the square fire source field, by comparing the results of experiments in 1:2 and full-scale tunnels. Hu et al. [39] studied, experimentally and computationally, the maximum temperature on the tunnel roof. They first compared the FDS results with the experimental measurements, and then, with the values, they calculated with the model of Kurioka, Oka, Satoh, and Sugawa [38]. The temperature predicted by the

FDS was found to converge quite well with both the experimental measurements and the results of the empirical model.

The Li et al. [40] study provided a theoretical analysis of the maximum temperature in the tunnel's roof, based on the axon-symmetric fire plume theory, taking into account the heat release rate, the longitudinal ventilation velocity, and the geometry of the tunnel. They conducted scale tests to investigate the maximum temperature on the tunnel roof in case of fire. Li and Ingason [41] analyzed the effects of different ventilation systems, ventilation speeds, heat release rates, tunnel geometries, and fire sources on the maximum temperature under the tunnel roof for large fires. They used data from many model-scale tests and largescale tunnel fire experiments around the world. They proposed some empirical formulas for the maximum temperature on the tunnel's roof nearby the fire source for low and high ventilation speeds. Rosignuolo et al. [42] merged the empirical relationships from the two studies above and defined a formula that exhibited non-linear regression and calculated temperature values as a function of distance from the fire source. They also distinguished between small and large fires because, in the event of a very large fire, the flame reaches the roof, so that the maximum temperature corresponds to the temperature of the flame, and it is difficult to estimate the mass flow rate. If a portion of the flame volume containing the combustion zone reaches the roof of the tunnel, the maximum temperature on the roof is reaching a constant value. Therefore, the maximum roof gas temperature is limited to a maximum of 1350 ◦C. The formula that calculates the temperature as a function of distance from the fire source is expressed as follows:

$$\frac{\Delta T(x)}{\Delta T\_{\text{max}}} = 0.55 \, \exp\left(-0.143 \frac{\chi - \chi\_{\upsilon}}{H}\right) + 0.45 \, \exp\left(-0.024 \frac{\chi - \chi\_{\upsilon}}{H}\right) \tag{9}$$

where *χ* is the distance from the fire source, *H* is the height of the tunnel, and *χν* is the virtual origin, which is calculated as

$$\chi\_{\vee} = \begin{cases} L\_f - 10H, & L\_f > 10H \\ 0, & L\_f \le 10H \end{cases} \tag{10}$$

where *Lf* is the flame length (m). The maximum temperature of the gas depends on its maximum value on the roof (Δ*Tmax*), which in turn was calculated by defining two different areas, depending on the dimensionless ventilation speed *V* :

$$
\Delta T\_{\text{max}} = \begin{cases}
17.5 \frac{\dot{Q}^{2/3}}{\dot{H}\_{ef}^{5/3}}, & V \le 0.19 \\
\frac{\dot{Q}}{\dot{u}\_0 b\_{f0}^{1/3} \dot{H}\_{ef}^{5/3}}, & V' > 0.19
\end{cases} \tag{11}
$$

where . *Q* is the total heat release rate (kW), *He f* is the effective height of the tunnel (the height from the fire source to the ceiling) (m), *u*<sup>0</sup> is the ventilation velocity, *bf* <sup>0</sup> is the radius of the fire source (m), and *V* is defined as:

$$V' = \frac{u\_0}{w^\*} \tag{12}$$

where *w*∗ is the characteristic velocity of the plume and is expressed by the following equation:

$$w^\* = \left(\frac{g\dot{Q}}{b\_{f0}\rho\_0 \mathbb{C}\_p T\_0}\right)^{1/3} \tag{13}$$

where *g* is the gravitational acceleration, (m/s2), *bf* <sup>0</sup> is the radius of the fire source (m), *ρ*<sup>0</sup> is the density of air (kg/m3), *Cp* is the specific heat capacity (kJ/kgK), and *T*<sup>0</sup> is the ambient temperature (K).

The values calculated with the empirical model were compared to the numerical results for Scenario C, measured with thermocouples at the tunnel ceiling (z = 7 m) above and downwind the fire source (where x = 0 is directly above the fire source). A large difference in temperature values is shown in Figure 4. The figure shows the good agreement of the curves at 0 m during the initial stages of the fire. From this point on, the temperature decreases and follows a scattering trend. There is a consistent difference of around 400 ◦C between the empirical and the actual temperature values for every location, except for above the fire pool. This difference leads to the assumption that, since the flame length increases up to the point that it is equal to and greater than the effective height of the tunnel and starts crawling horizontally to the ceiling, the smoke layer surrounds the fire source. The results indicate that there was not enough oxygen for the fire to grow further, so the fire became ventilation controlled. Other possible reasons for these differences might be that there are limitations to the validity of parameters such as the geometry of the tunnel and the fire area (the FDS fire source has different dimensions than those used in the experiments).

**Figure 4.** Temperatures at the tunnel ceiling (z = 7 m) above (x = 0 m) and downwind the fire source (x = 10, 30, 50 and 100 m) for Scenario C, calculated using the empirical model and the actual temperatures from the FDS output.

Additionally, for the validation of temperature, the results were compared with the hydrocarbon (HC) fire curve and modified hydrocarbon fire curve (HCM) used for fire resistance testing in structural designs and the RABT-ZTV train fire curve. The HC fire curve is used when there is a risk of liquid fires, e.g., in industrial plants (chemicals), on ships, pontoons, oil platforms, road, and railway tunnels, as well as for petrol and diesel pool fires. The HCM fire curve is derived from the HC curve and adjusted for French regulations. Both HC and HCM curves showed an abrupt temperature gradient in the first minutes, and their difference is the maximum temperature of 1100 and 1300 ◦C, respectively. The guideline for the equipment and operation of road tunnels (Richtlinie

für die Ausstattung und den Betrieb von Straßentunneln) (RABT) defines the so-called RABT tunnel fire. This is often used for the fire protection design of tunnel structures. The RABT tunnel fire curve is an artificial fire curve and is intended to cover the maximum temperatures resulting from tunnel fires [19]. The curve can be used in tunnel design and testing. The HC and HCM fire curve equations are shown below in Equation (14) and Equation (15), respectively.

$$T = 20 + 1080 \cdot \left(1 - 0.325 \cdot e^{-0.167 \cdot t} - 0.675 \cdot e^{-2.5 \cdot t} \right) \tag{14}$$

$$T = 20 + 1280 \cdot \left(1 - 0.325 \cdot e^{-0.167 \cdot t} - 0.675 \cdot e^{-2.5 \cdot t} \right) \tag{15}$$

where *t* is time in minutes.

The RABT-ZTV train fire curve followed a linear growth for the first 5 min, where it reached a maximum temperature of 1200 ◦C. A comparison of the ceiling gas temperature from the FDS output for Scenario C and the HC, HCM, and RABT-ZTV curves is shown in Figure 5. The FDS output was in good agreement with the RABT-ZTV fire curve for the first minutes, until the activation of the ventilation system, which happened at 480 s. Scenario C was selected because of the delay in the activation of the jet fans, so the temperature rise was more in accordance with the fire curves.

**Figure 5.** Tunnel ceiling temperature (z = 7 m) above the fire source (x = 795 m), and comparison of the FDS output for Scenario C and the HC, HCM, and RABT ZTV standard fire curves.

### *2.3. Evacuation and Tenability Analysis*

For each scenario (A, B, and C) three different evacuation times were examined +60 s, +120 s, and +300 s after the start of the fire. The evacuation of passengers was effectuated along the walking platform on the side of the tunnel. The total number of passengers was 524 which corresponds to 80% capacity of an intercity train with 8 wagons, each of which can accommodate 80 passengers. The evacuation simulation and FED index calculation were conducted with Pathfinder, which is an interface for the FDS + EVAC code developed by the NIST [31].

For the three evacuation scenarios, eight occupant sources were created, each one simulating the door of a wagon. The total population number that exited the tunnel was 524 persons. This accounts for approximately 80% capacity of the train, assuming each wagon carries 80 passengers at full capacity. The profile of the passengers was set as the default profile provided, with a height of 1.8 m and a walking speed of 1.19 m/s. For the sake of convenience, only the entrance of the tunnel was marked as the exit point for the passengers to move towards to.

#### **3. FDS Simulation Results**

Three different scenarios were examined (A, B, and C) with three different fan activation times after a fire accident, at +180 s, +300 s, and +480 s, respectively. For all the cases, all the fans were activated at the same time. The three different scenarios for the different jet fan activation times were set up with the Pyrosim interface provided by Thunderhead Engineering. The extracted results were temperature, flow velocity, the volume fraction of gases, and heat release rate (HRR) values.

As shown in Figure 6, the input fire curve is in good agreement with the numerical results for the HRR curve.

Figures 7 and 8 show that the longer the activation of the ventilation system is delayed, the higher the temperatures appear in the tunnel ceiling. This rapid rise in high temperature can result in spalling, which is the detachment of structural material, such as concrete, often happening in an explosive way. This poses a great threat for the population in the tunnel.

Moreover, a time delay in the ventilation system activation creates a greater smoke dispersion, as shown in Figure 9, and thus a greater dispersion of the harmful combustion products. This can affect the evacuation process due to reduced visibility and soot inhalation.

As shown in Figure 10, the temperatures at the ceiling above the fire source are in good agreement with the RABT ZTV fire curve in the first minutes of the simulation. Then, when the ventilation system was activated, the temperature dropped and reached a steady value of around 450 ◦C.

Figure 11 shows that the development of ventilation flow started immediately after the activation of the jet fans. The air flow inside the tunnel needs about 300 s, after the jet fans achieve the nominal power, to reach the critical velocity value, but the thrust of the jet fans, even before they reach full power, is capable of pushing the smoke downwind of the fire source. Moreover, Scenario A showed some smoke layer entertainment, which was due to the pressure difference at the tunnel openings, which created a draft with a velocity around 0.5 m/s (Figure 12 at the 180 s time frame). As shown in Figure 12, the peak at 750 m (fire source) was due to the smoke layer and the presence of the train, which reduced the area of the tunnel in this section, resulting in higher velocities. This phenomenon can be seen in Figure 11, where, near the train location, a great volume of the flow had the critical velocity value. As depicted in Figure 13, the mean longitudinal velocity was in the range of 2.7 m/s and was achieved within 600 s from the jet fan activation time. Ventilation conditions are achieved as long as the jet fans reach their nominal output volume flow, but this does not mean the tunnel airflow velocity has reached the critical value.

**Figure 6.** FDS input fire curve and FDS results fire curve.

**Figure 7.** Tunnel ceiling (z = 7 m) gas temperature above (x = 0 m) and downwind (x = 10, 30, 50, and 100 m) from the fire source (**a**) for Scenario A, (**b**) Scenario B, and (**c**) Scenario C.

**Figure 8.** Tunnel ceiling (z = 7 m) temperatures for each scenario at the jet fan activation time frame: (**a**) Scenario A at +180 s, (**b**) Scenario B at +300 s, and (**c**) Scenario C at +480 s after the fire accident.

**Figure 9.** Smoke dispersion along the tunnel and backlayering at the jet fan activation time frame for (**a**) Scenario A at +180 s, (**b**) Scenario B at +300 s, and (**c**) Scenario C at +480 s.

**Figure 10.** FDS ceiling temperatures above the fire source (x = 795 m, z = 7 m) for all three scenarios in comparison with the RABT ZTV fire curve.

**Figure 11.** Development of the critical ventilation velocity of 2.7 m/s (isosurface with the green colour), at the train location, for (**a**) +90 s, (**b**) +180 s, and (**c**) +300 s after the activation of the jet fans. The jet fans reached their maximum output 90 s after their activation, while the development of the flow was reached approximately after 300 s.

**Figure 12.** Lengthwise U velocity at symmetry plane y = 5.5 m and z = 3.5 m for Scenario A at different time frames (t = +180 s, +270 s, +600 s, and +720 s).

**Figure 13.** Lengthwise mean U velocity for Scenario A, upwind from the fire source.

### **4. Temperature and FED Index**

Increased air temperatures may make it difficult for the people to move. This increase affects the train evacuation process. The critical temperature in a confined space, such as a tunnel, is the temperature where a person's movement and actions are impaired, and it is defined as 60 ◦C [43]. Figure 15 shows that, for every scenario, the temperature at the walking platform, at a height of 2 m, did not exceed 24 ◦C. In these conditions, the evacuation can be carried out normally.

The calculation of FED with Pathfinder used the equations described in the SFPE Handbook of Fire Protection Engineering, 5th Edition [44]. Only the concentrations of the narcotic gases CO, CO2, and O2 were used to calculate the FED value.

$$\text{FED} = \text{FED}\_{\text{CO}} \cdot \text{V}\_{\text{CO}\_2} + \text{FED}\_{\text{O}\_2} \tag{16}$$

This calculation does not include the effect of hydrogen cyanide (HCN), and the effect of CO2 is only due to hyperventilation. Carbon dioxide does not have toxic effects at concentrations of up to 5%, but it stimulates breathing, which increases the rate at which the other fire products are taken up.

$$\text{FED}\_{\text{CO}} = 3.317 \cdot 10^{-5} \cdot \text{CO}^{1.036} \cdot V \cdot t/D \tag{17}$$

where CO is the concentration of carbon monoxide (ppm *v*/*v* 20 ◦C), *V* is the volume of air breathed each minute (L/min), with its value being 25 L/min for an activity level of light work (walking to escape), *t* is time in minutes, and *D* is the exposure dose (% COHb), with its value being 30% for an activity level of light work.

Hyperventilation due to carbon dioxide can increase the rate at which fire products are taken up. The multiplication factor is given by the equation

$$V\_{\rm CO\_2} = \exp(0.1903 \cdot \% \text{CO\_2} + 2.0004)/7.1\tag{18}$$

where % CO2 is the volume fraction of CO2 (*v*/*v*). The fraction of an incapacitating dose of low O2 hypoxia is calculated as:

$$\text{FED}\_{\text{O}\_2} = \frac{t}{\exp\left[8.13 - 0.54 \cdot (20.9 - \%\text{O}\_2)\right]}\tag{19}$$

where *t* is time in minutes, and % O2 is the oxygen volume fraction (*v*/*v*). The measurement location for FED quantity sampling is 90% of occupant height. For example, FED data for an occupant with height = 1.8 m would be sampled 1.62 m above the occupant's location. When the occupant is not inside an FDS mesh, the FED calculation is paused until the occupant enters another mesh [31].

FED values never reached the critical value of 1, where a person is considered incapacitated, but some interesting results were extracted. The maximum values of FED are presented in Table 1 and compared in Figure 14. As shown, the longer the activation of the ventilation system was delayed, the more the uptake of combustion products and toxic gases increases. This is also shown in Table 2 and Figure 16, where the results for the mean fed values for each case are presented. It is important to be mentioned that, while maximum FED values for scenario A are about six times lower than the values for Scenario C, the mean FED values for Scenario A are an order of magnitude lower than Scenario C.

**Figure 14.** Max FED index for the three different scenarios and evacuation times.


**Table 1.** Max FED index for the three different scenarios and evacuation times.

**Table 2.** Mean FED index for the three different scenarios and evacuation times.


**Figure 15.** Temperature distribution upwind from the fire source at y = 0.5 m (middle of the walking platform) and at a height of 2 m for (**a**) Scenario A, (**b**) Scenario B, and (**c**) Scenario C. Temperatures do not exceed 24 ◦C.

**Figure 16.** Mean FED index for the three different scenarios and evacuation times.

### **5. Discussion and Conclusions**

Three different scenarios with different jet fan activation times, combined with three different evacuation scenarios, were studied. Firstly, validation of the temperature and critical velocity was conducted. The results showed a good agreement between the theoretical and simulation values, thereby the ventilation system was capable of providing sufficient air for smoke extraction.

The evacuation process showed that the passengers were not in danger, but important results regarding the importance of the jet fan activation time were extracted: From the FED values calculated for the nine cases in total, it was concluded that the most important action in a tunnel fire is the ventilation system activation time, after the start of a fire. If the ventilation system is delayed for any reason and the passengers have started evacuating the train, it is safer to wait for the population to travel some distance, possibly way past the train, and then the jet fans should be activated. If the ventilation system is activated in time, it seems safer to delay the initialization of the evacuation until the jet fans have reached their nominal output.

It is important to mention that there was no diversity in the population. Each simulated person had the same characteristics with each other as the movement speed and other factors that affect the evacuation process such as dimensions and comfort distance from each other. This affects results since, in real-life scenarios, the population is mixed (elderly people, people requiring assistance to move, people with baggage and families).

**Author Contributions:** Conceptualization, T.Z. and K.V.; methodology, T.Z., K.V. and I.S.; software, T.Z. and K.V.; validation, T.Z. and K.V.; formal analysis; T.Z.; resources, K.V. and I.S.; data curation, T.Z.; writing—original draft preparation, T.Z. and K.V.; writing—review and editing, K.V. and I.S.; visualization, T.Z.; supervision, K.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used in this study are available upon request from the corresponding author. The data can be easily reproduced from the theoretical analysis described in the study.

**Acknowledgments:** This work was supported by computational time granted from the National Infrastructures for Research and Technology S.A. (GRNET S.A.) in the National HPC facility—ARIS under project ID pr011045-UrbanFirePlan2. The authors acknowledge and express particular appreciation to Thunderhead Engineering for providing licenses for Pyrosim and Pathfinder software, which was used for the designing, simulation scenarios and visualization of the results.

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

### **References**


### *Article* **Air Flow Study Around Isolated Cubical Building in the City of Athens under Various Climate Conditions**

**Chariton L. Pavlidis, Anargyros V. Palampigik, Konstantinos Vasilopoulos, Ioannis C. Lekakis and Ioannis E. Sarris \***

> Laboratory of Thermo Fluid Systems (LTFS), Department of Mechanical Engineering, Ancient Olive Grove Campus, University of West Attica, Thivon Str. 250, Egaleo, 12244 Athens, Greece; cpavlidis@uniwa.gr (C.L.P.); msrtf21X03@uniwa.gr (A.V.P.); kvassil@uniwa.gr (K.V.); lekakis@uniwa.gr (I.C.L.) **\*** Correspondence: sarris@uniwa.gr; Tel.: +30-2105381131

**Abstract:** This study focuses on the airflow and pollutant dispersion around an isolated cubical building located in a warm Mediterranean climate, taking into account the local microclimate conditions (of airflow, albedo of building and soil, and air humidity) using a large-eddy simulation (LES) numerical approach. To test the reliability of computations, comparisons are made against the SILSOE cube experimental data. Three different scenarios are examined: (a) Scenario A with adiabatic walls, (b) Scenario B with the same constant temperature on all the surfaces of the building, and (c) Scenario C using convective and radiative conditions imposed by the local microclimate. For the first two cases the velocity and temperature fields resulting are almost identical. In the third case, the resulting temperature on the surfaces of the building is increased by 19.5%, the center (eye) of the wake zone is raised from the ground and the maximum pollutant concentration is drastically reduced (89%).

**Keywords:** LES; microclimate model; thermal radiation; pollutant dispersion; urban planning

### **1. Introduction**

The excessive population concentration in megacities has resulted in environmental pollution problems and overcrowded living conditions. The air quality in these cities is determined by the airflow in the complex urban terrain, the temporal and spatial conditions of pollutant sources, and the local urban microclimate parameters. The airflow and pollutant dispersion in an urban environment can be studied in different geometrical scales, such as blocks of buildings, street canyons, and isolated buildings. Several wind-tunnel, real-scale experiments and numerical studies examine the airflow and pollutant dispersion around blocks of buildings [1–5] and street canyons [6–12].

The city's local microclimate parameters are air and surface temperatures, humidity, direct and diffuse solar irradiation, and wind speed at different directions [13]. The knowledge of the urban microclimate is necessary to prevent and reduce pollution problems. The local microclimate conditions are playing an important role in the airflow and pollution distribution in an urban environment [14]. For this reason, complex models can be used to calculate mandatory variables such as the heat capacity of a building, the shortwave and longwave radiation according to the solar position, the emissivity, and the air humidity. According to these variables, the increment or decrement of the thermal comfort rate can be found to define the human comfort conditions. Another important issue is the values of the heating rate and of the thermal absorption of the surfaces of the building. Using complicated microclimate models, the distribution of the environmental temperature and their effect on pedestrians can be defined [15–19]. Small-scale meteorological models are widely used to define the comfort conditions inside the urban environment for winter and summer seasons [20–23].

**Citation:** Pavlidis, C.L.; Palampigik, A.V.; Vasilopoulos, K.; Lekakis, I.C.; Sarris, I.E. Air Flow Study Around Isolated Cubical Building in the City of Athens under Various Climate Conditions. *Appl. Sci.* **2022**, *12*, 3410. https://doi.org/10.3390/ app12073410

Academic Editor: Rüdiger Schwarze

Received: 18 January 2022 Accepted: 21 March 2022 Published: 27 March 2022

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

**Copyright:** © 2022 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/).

A phenomenon by which warm layers of air are accumulated in densely populated centers is defined as an Urban Heat Island (UHI) creating uncomfortable thermal conditions for urban living [24]. A high air-temperature increment leads to the need for higher energy consumption for cooling in the summer, with a consequence of increasing air pollutants in the atmosphere [25].

It is also identified that the rapid augmentation of the air temperature depends on the airflow velocity. At calm wind conditions, where the airflow velocities are lower than 3 m/s, the density of hot gas masses is much higher than at heavy wind gusts [26]. As a result, wind gusts with high airflow velocities can remove the air hot masses from the urban environments, which leads to the decrement of the UHI phenomenon [27].

The effect of the microclimate in densely populated urban centers is also an important topic for study. Gronemeier and Raasch [28] studied the urban flow for the city of Hong Kong with PALM software. The numerical results showed a rapid increment in wind velocity at the pedestrian crossing areas and high changes of the airflow direction [28,29]. Another study by Pfafferott and Rißmann [30] with the PALM code examined the interaction between buildings and the urban microclimate using an energy balance solver. They found that the temperature increment of the urban environment was caused by the heated building, as the outcome of energy from the indoor model [30]. Resler and Eben [31] examined the airflow over the city of Prague in the Czech Republic with PALM code. Their study showed that the temperature and the concentration of NOx are in a good agreement with real field measurements during the summer season. In the winter season, the deviations between numerical results and experimental data pointed to inaccuracies in modeling the atmospheric boundary layer [31]. Thus, real field urban measurements can give important information for both real-scale meteorological phenomena and their impact on an urban boundary layer. To better understand the complex phenomena of the airflow in an urban environment, it is important to combine real field measurements and wind tunnel experiments so that these two methods can be interrelated [32,33]. Numerical methods can be also used to simulate the turbulence modelling of the airflow around both small- and large-scale cubical obstacles [34–38].

Several wind-tunnel experiments exist that describe the airflow characteristics around a cubical building. The wind-tunnel experiments are suitable for collecting detailed spatial experimental data [32,33,39–41] under stable conditions. However, experimental wind tunnel experiments are not efficient for the prediction of real-scale meteorological conditions.

The present study is focusing on the simulation of the flow and pollution distribution around an isolated building [42] and, more importantly, on the influence of a meteorological model on these results, i.e., what differences exist between the use of real boundary conditions and those employed for reasons of convenience. In addition, the use of an isolated building instead of a large urban area ensures good reliability of the numerical results for the available resources. In this respect, such studies do not seem to exist in the scientific literature.

The isolated building of this study is in the climatological area of Athens, Greece. Three different scenarios are examined. Scenario A studies the neutral flow and pollutant dispersion in the wake area of an isolated cubical building with adiabatic walls, without any buoyancy forces. Scenario B studies the same building with specified constant temperature of 50 ◦C on its surfaces. Scenario C studies the same building with convective and radiative conditions imposed by the local microclimate.

### **2. Configuration and Pollutant Dispersion Modelling**

### *2.1. Problem Configuration*

Figure 1 presents the geometry of the computational domain around the isolated cubical building of this study, with exactly the same dimensions as the SILSOE cube. The SILSOE cube field experiment for the study of the airflow around it [43], with a 6 m height cube placed in an open area, has been used as a test case in several other studies [44]. This geometry remains unchanged in all the different simulated scenarios of this work.

**Figure 1.** Computational domain and boundary conditions.

The global coordinate system origin is the frontal left edge of the computational domain and the distance between the frontal area of the computational domain and the building is 5 H. The distance between the lateral boundaries of the computational domain and the building's surface is also 5 H [45]. The distance from the rear surface of the building up to the outlet boundary of the computational domain is set at 12 H and the total height of the computational domain is 5 H. The resulting blockage effect value is approximately 1.81%. Under the recommendations of the German Association of Engineers (VDI), the blockage effect should be maintained below 10% during the simulations [46] and according to [46,47] should be even lower than 3%, which is satisfied in the present study.

The Reynolds number is kept constant at 4.03 × 106, based on the height of the cube and the free stream velocity. The specified flow field at the inlet of the computational domain is in the form of a logarithmic profile with the no-slip condition at the ground. A passive source of pollutants is located on the floor of the computational domain at a distance *H* from the rear surface of the building, using as a pollutant methane (CH4) gas with a concentration and release of 600 ppm and 18.5 L/h, respectively. The dimensionless concentration coefficient *K* of the passive pollutant is defined as [48,49]:

$$K = \frac{(\mathbf{c}\_{Measured} / \mathbf{c}\_{Source})\mu\_H H^2}{Q\_{Source}} \tag{1}$$

where *cMeasured* is the tracer concentration at any position, *cSource* is the tracer concentration on the pollutant source, *QSource* is the release rate of the pollutant, and *uH* is the velocity of the flow at the height of the building.

In Scenarios B and C, the temperature is affecting the flow field through the usual Boussinesq approximation. In these cases, the Grashof number, i.e., the ratio of thermal buoyancy forces over viscous forces is defined as [50]:

$$Gr\_H = \frac{g\beta \left(T\_{Surface} - T\_{\infty}\right)H^3}{v^2} \tag{2}$$

where *TSur f ace* is the temperature on the surfaces of the building, *T*<sup>∞</sup> is the bulk temperature, *β* is the coefficient of thermal expansion, and *v* is the kinematic viscosity of the fluid. In Scenarios B and C, the Grashof number based on the height of the building is 2.43 × 1011 and 7.68 × 1011, respectively. Thus, the buoyancy forces are intensive and the boundary layers close to the heated surfaces are highly turbulent.

#### *2.2. Governing Equations*

The PALM model system, an open-source software that simulates atmospheric and oceanic boundary layers, is used for the present simulations. The Navier-Stokes equations in a non-hydrostatic, filtered, and incompressible form under the Boussinesq approximation are solved. The basic equations for the conservation of mass, momentum, energy, and moisture are defined as:

$$\frac{\partial \mu\_i}{\partial x\_j} = 0 \tag{3}$$

$$\frac{\partial u\_i}{\partial t} = -\frac{\partial u\_i u\_j}{\partial x\_j} - \varepsilon\_{ij\bar{k}} f\_j u\_k + \varepsilon\_{i\bar{3}j} f\_3 u\_{\bar{\varrho},j} - \frac{1}{\rho\_0} \frac{\partial \pi^\*}{\partial x\_i} + g \frac{\theta\_v - \langle \theta\_v \rangle}{\langle \theta\_v \rangle} \delta\_{l3} - \frac{\partial}{\partial x\_j} \left( \overline{u\_i'' u\_j''} - \frac{2}{3} \epsilon \delta\_{l\bar{j}} \right) \tag{4}$$

$$\frac{\partial \theta}{\partial t} = -\frac{\partial u\_j \theta}{\partial x\_j} - \frac{\partial}{\partial x\_j} \left( \overline{u\_j'' \theta''} \right) - \frac{L\_v}{\mathbb{C}\_p \Pi} \Psi\_{q\_\upsilon} \tag{5}$$

$$\frac{\partial q\_{\upsilon}}{\partial t} = -\frac{\partial \mu\_{j} q\_{\upsilon}}{\partial x\_{j}} - \frac{\partial}{\partial x\_{j}} \left( \overline{\mu\_{j}^{\prime\prime} q\_{\upsilon}^{\prime\prime}} \right) + \Psi\_{q\_{\upsilon}} \tag{6}$$

$$\frac{\partial \mathbf{s}}{\partial t} = -\frac{\partial \mathbf{u}\_j \mathbf{s}}{\partial \mathbf{x}\_j} - \frac{\partial}{\partial \mathbf{x}\_j} \left( \overline{\mathbf{u}\_j'' \mathbf{s}''} \right) + \Psi\_s \tag{7}$$

where *i*, *j*, *k* ∈ {1, 2, 3}, the components of the velocity (*u*<sup>1</sup> = *u*, *u*<sup>2</sup> = *v* and *u*<sup>3</sup> = *w*) are defined by the *ui* variable at specific positions in the flow field *xi*, where *x*<sup>1</sup> = *x*, *x*<sup>2</sup> = *y* and *x*<sup>3</sup> = *z*, and *t* is time. *qv* and *s* are the latent heat transfer and moisture, respectively. The gravitational acceleration is denoted by *g*, the density of the dry air is defined by *ρ*0, and *Lv* is the latent heat of vaporization. The source term of the variable *qv*, is defined as Ψ*qv* and the sink term of the variable *s* is defined as Ψ*s*. The angle brackets correspond to the horizontal averages of the flow field and the double prime indicated the subgrid-scale variables.

The subgrid-scale turbulence kinetic energy is defined as:

$$
\sigma = \frac{1}{2} \left( \overline{u\_i'' u\_i''} \right) \tag{8}
$$

The modified perturbation pressure based on the perturbation pressure *p*∗ is defined as:

$$
\pi^\* = p^\* + \frac{2}{3}\rho\_0 c \tag{9}
$$

The potential temperature is defined from the equation:

$$\theta = \frac{T}{\Pi} \tag{10}$$

where *T* is the instantaneous absolute temperature and the Exner function is calculated by the equation:

$$
\Pi = \left(\frac{p}{p\_0}\right)^{\frac{p\_d}{\ell\_p}} \tag{11}
$$

where *p* is the hydrostatic pressure of the air, *p*<sup>0</sup> = 1000 hPa is the reference pressure, *Rd* is the specific gas constant for dry air, and *Cp* is the specific heat capacity under constant pressure.

The virtual potential temperature is defined by the equation:

$$\theta\_{\upsilon} = \theta \left[ 1 + q\_{\upsilon} \left( \frac{R\_{\upsilon}}{R\_d} - 1 \right) - q\_l \right] \tag{12}$$

where *Rv* is the specific gas constant for water vapors, and *ql* is the liquid water mixing ratio.

### *2.3. Urban Surface Model (USM) of the PALM Model System*

The urban surface model (USM) of the PALM model system is used for the energy balance in the computational domain. The energy balance solver is driven by three main procedures: (a) the solver predicts the temperature of the outer surfaces, (b) the turbulence dispersion of the sensible temperature is calculated at the surfaces near the walls, and (c) the calculation of the sublayer heat flux caused by convection. The first two procedures of the energy balance solver are executed simultaneously. The third procedure is executed under the calculations of a subgrid-scale model that predicts the thermal diffusion from the bluff body.

The energy balance solver is defined as:

$$C\_0 \frac{dT\_0}{dt} = R\_n - H - LE - G \tag{13}$$

where *C*<sup>0</sup> is the specific heat capacity, *T*<sup>0</sup> is the radiative temperature of the surface skin layer, *Rn* is the net radiation, *H* is the sensible heat, *LE* is the latent heat, and *G* is the heat flux at the surface of the ground. On the USM, the energy balance is calculated separately for each surface, and the three different types of radiations (sensible heat, latent heat, and ground heat flux) from the surface heat are combined.

The calculation of the convective heat transfer between the air and the outer surfaces is defined by the equation:

$$H = h(\theta\_1 - \theta\_0) \tag{14}$$

where *θ*<sup>0</sup> is the temperature of the outer surface of the building, *θ*<sup>1</sup> is the temperature of the air mass in contact with the outer surface of the building, and *h* is the coefficient of convective heat transfer and is parameterized only for the vertical surfaces inside the computational domain [51]. For the horizontal surfaces inside the computational domain, the coefficient of convective heat transfer is parameterized under the Monin-Obukhov similarity, which involves the calculation of the local friction velocity [52,53]. The friction velocity is used for the calculation of flow momentum near the surface, for each surface separately.

### **3. Initial and Boundary Conditions**

The inlet velocity profile is described by the logarithmic law of the wall with a freestream velocity value of 10.13 m/s. Figure 2 presents the normalized velocity profiles at both the inlet of the computational domain and in the X/H = 2 position from the inlet for Scenarios A–C, where A and B should be identical.

**Figure 2.** Normalized velocity profiles at the (**a**) inlet of the computational domain and (**b**) X/H = 2 position from the inlet for Scenarios A−C [44,54].

As observed in Figure 2a, an important deviation appears in Scenario C that corresponds to the application of the meteorological model in comparison with the other two scenarios. This difference appears also in the free-stream flow region and is caused by the strong buoyancy effects caused by the increased thermal heating of the air.

As shown in Figure 3, the present numerical results for the pressure coefficient are in good agreement with the corresponding experimental data of Scenarios A and B.

**Figure 3.** Pressure coefficient distribution at the building's symmetry plane for Scenarios A−C [4,44,54].

The pressure coefficient at the surfaces of the building for all scenarios is defined as:

$$C\_P = \frac{p - p\_{REF}}{\frac{1}{2}\rho\_{REF}\mu\_{\infty}^2} \tag{15}$$

where *p* is the static pressure of the fluid at a point, *pREF* and *ρREF* are the static pressure and density of the fluid at the free-stream at the inlet of the computational domain, respectively, and *u*∞ is the free-stream velocity at the inlet.

As observed in Figure 3, important differences exist between Scenarios A–C because of the heating of the surfaces of the building and the microclimate model application. The pressure coefficient distribution on the frontal surface of the building for Scenarios B and C is higher than the corresponding for Scenario A because of the influence of the mixed convection heat transfer, which causes reduction on the density of the streamlines in this region. In Scenario C, this influence is expected to be larger than for Scenario B. In Scenario B, the pressure distribution at the rear surface of the building significantly approaches the corresponding for Scenario A. In Scenario C, the pressure distribution at the rear surface of the building is the highest compared with other scenarios because of the domination of the buoyant heating forces there.

Figure 4 presents the vertical profile of the turbulent kinetic energy at the inlet of the computational domain for Scenarios A–C, as compared with the corresponding numerical data on the same normalized upstream position from the cube [18].

As described in Figure 4, the turbulent kinetic energy distribution at the inlet of the computational domain is almost identical for Scenarios A and B and different from Scenario C because of the effect of the buoyancy forces from the microclimate conditions.

**Figure 4.** Turbulent kinetic energy at the inlet of the computational domain for Scenarios A–C [18,44].

On the lateral boundaries of the computational domain, periodic boundary conditions are employed because of the temporal and spatial periodicity of the flow. Neumann boundary conditions are used for the turbulent kinetic energy (*e*), temperature (*θ*), and the perturbation pressure (*p*∗) to calculate the following equations concerning the building's height:

$$
\overline{\varepsilon}\left(-\frac{\Delta y}{2}\right) = \overline{\varepsilon}\left(+\frac{\Delta y}{2}\right) \tag{16}
$$

$$
\overline{\theta}\left(-\frac{\Delta y}{2}\right) = \overline{\theta}\left(+\frac{\Delta y}{2}\right) \tag{17}
$$

$$
\overline{p^\*}\left(-\frac{\Delta y}{2}\right) = \overline{p^\*}\left(+\frac{\Delta y}{2}\right) \tag{18}
$$

The Dirichlet boundary condition is used for the velocity at the ground of the computational domain for the no-slip condition. Therefore, the velocity components at the ground are:

$$
\overline{u}(z=0) = 0\tag{19}
$$

$$
\overline{v}(z=0) = 0\tag{20}
$$

$$
\overline{w}(z=0) = 0\tag{21}
$$

On the staggered computational grid, the velocity components in the X and Y directions, respectively, are defined at the specific height of *<sup>z</sup>* <sup>=</sup> <sup>±</sup>Δ*<sup>z</sup>* <sup>2</sup> . Thus, the symmetry boundary condition is used as:

$$
\pi \left( -\frac{\Delta z}{2} \right) = -\pi \left( +\frac{\Delta z}{2} \right) \tag{22}
$$

$$
\overline{\upsilon}\left(-\frac{\Delta z}{2}\right) = -\overline{\upsilon}\left(+\frac{\Delta z}{2}\right) \tag{23}
$$

On the upper surface of the computational domain, the perturbation pressure (*p*∗) obeys the Neumann boundary condition. For the velocity to maintain the free-stream regime over the atmospheric boundary layer, Neumann boundary conditions are also applied to the velocity components of the velocity in the X and Y directions:

$$
\left.\partial\_z \overline{u}\right|\_{top} = \text{constant} \tag{24}
$$

$$\left. \partial\_{\overline{z}} \overline{v} \right|\_{top} = \text{constant} \tag{25}$$

Additionally, on the upper surface of the computational domain, Dirichlet boundary conditions are used for both the turbulence kinetic energy (*e*) and the temperature (*θ*).

#### **4. Nested Computational Grid**

PALM provides the capability of computational grid self-nesting. In the self-nesting mode, there is a root/parent computational domain with the ability of nesting up to 63 levels of nested/child computational subdomains. Each subdomain may be root/parent to another nested subdomain. As a result, the second subdomain is simultaneously nested referring to the first parent, and also root, to the third nested. The nested computational subdomain receives all the appropriate and mandatory information for the three components of the velocity and all the prognostic variables of the vectors from its boundaries with the root computational domain. The flow field data are interpolated from the coarse to the finer computational grid. By the end of each time step, the corrected solution from the solver is reversely interpolated to the root/parent domain.

Figure 5 shows the grid arrangement used in the simulations. The root computational grid (coarse grid) is shown at full height of the computational domain in contrast with the nested (finer grid-darker parts), which is shown at two building heights.

**Figure 5.** Root and nested computational domains of the flow field.

For the simulations, three different computational grids are used with an incremental resolution from the coarse to the finer: the coarse computational grid with 924,400 cells, the medium grid with 1,248,400 cells, and the fine grid with 1,842,400 cells. The grid errors using two computational grids are estimated from the following equation [4]:

$$GCI = \frac{f\_2 - f\_1}{1 - r^p} \tag{26}$$

where *f*<sup>2</sup> is the numerical solution that results from the medium computational grid and the corresponding from the finer is *f*1. *r* is the refinement factor between the two computational grids and *p* is the accuracy of the algorithm, which is 3 for the present study. The refinement factor between the medium and finer computational grid is approximately 1.47.

Figure 6a presents the profile of the error bars of the normalized velocity, at the axial position 0.5 H from the rear surface of the building for scenario A. Additionally, Figure 6b shows the profile of the normalized concentration of the pollutant at the same axial position and scenario.

The mean value of the errors of the normalized velocity between the coarse and the medium computational grids is approximately 3.77% and the corresponding value between the medium and the fine grids is about 1.54%. The mean value of the errors of the pollutant between the coarse and the medium computational grids is approximately 2.77% and the corresponding value between the medium and the fine grids is about 1.24%. For both of the aforementioned variables, the decrement of the errors going from the coarse to the finer computational grids is clearly observed.

**Figure 6.** GCI error bars estimated from the fine to medium grids for Scenario A: (**a**) for the normalized velocities ux/u∞ and (**b**) for the normalized concentration of the pollutant at the axial position of 0.5 H from the rear surface of the building.

### **5. Numerical Details**

The global discretization of the computational grid is accomplished using finite differences on a staggered Cartesian Arakawa-C grid [55]. The equations are spatially discretized on a fifth-order differential upwind scheme [56], while a third-order Runge-Kutta scheme is used on the temporal discretization [57]. The equations are implicitly filtered by the discretization of the computational grid and the subgrid-scale processes are calculated by the 1.5-order Deardorff numerical scheme [58]. Thus, it is assumed that the energy that is transported by the subgrid-scale vortices is proportional to the local gradients of the mean quantities [59,60].

The convergence criterion maintained on each simulation is kept below the 10−<sup>4</sup> value for each computed variable based on the error. The time step is automatically adjusted for the CFL constant value of 0.9.

Every simulation is carried out until flow stationarity has been achieved. Figure 7 presents the normalized velocity components with respect to the free-stream velocity (ux/u∞, uy/u∞, uz/u∞), at the point with Cartesian coordinates of X: 10 H, Y: 5.5 H, and Z: 0.5 H and for the period from 600 to 1400 s for Scenario A.

**Figure 7.** Velocity fluctuations of ux/u∞, uy/u∞, uz/u∞ at point X: 10 H, Y: 5.5 H, Z: 0.5 H from 600 up to 1400 (s) for Scenario A.

As it is shown in Figure 7, the flow seems to have attained stationarity conditions from the time point of 700 (s) onward. The attainment of stationarity conditions is confirmed by finding out that the mean value and higher order statistics are independent of the time of initiation of the measurements. If this is valid for the mean value and the autocorrelation function, the process is said to be weakly stationary. In this respect, the mean, the variance, and the autocovariance with a time delay of 2 (s) values are computed and given in tabulated form below (Table 1):

**Table 1.** Mean Value, Variance and Autocovariance for Specific Time Periods.


These computed values show that the flow field is at least weakly stationary.

However, although the flow field has reached stationary conditions it is possible the concentration field has not reached it yet.

Figure 8 presents the pollutant dispersion inside the computational domain for Scenarios A and C at the position with coordinates of X: 7.5 H, Y: 5.5 H, and Z: 0.5 H for the period of 200–1000 (s).

**Figure 8.** Concentration of the pollutant for Scenarios A and C at the position with coordinates of X: 7.5 H, Y: 5.5 H and Z: 0.5 H for the time period of 200 to 1000 (s).

The mean, the variance, and the autocovariance with time delay of 2 (s) values of the pollutant concentration are computed and given in tabulated form below (Table 2), which demonstrates that a passive scalar attains stationarity earlier in time than the velocity field:



### **6. Results and Discussion**

The air flow field around the cubical geometry is presented in Figure 9 using normalized velocity profiles on the symmetry plane at different streamwise positions upstream, downstream, and at the middle of the roof of the cube (X/H = 5.5) for all three scenarios [44,54,61–65]. Comparisons are also made with Richards and Castro's experimental data [60] for Scenario A.

**Figure 9.** Normalized velocity profiles on the symmetry plane at positions X/H = 0, X/H = 2, X/H = 4, X/H = 5.5, X/H = 7, X/H = 7.5, X/H = 8 and X/H = 10 for Scenarios A–C [44,60].

Differences in the normalized velocity profile between Scenario C and Scenarios A and B are observed especially in the wake region of the cube, where the mixing layer for Scenario C has moved almost one cube height above the mixing layer of the other two scenarios. This is caused by the strong buoyant forces created in this case.

Figure 10 shows the streamlines of the flow for Scenario A. A detachment of the flow at the top surface of the building and a recirculation of the flow on both the frontal and the rear surfaces of the building are observed.

Different recirculation lengths **Xf**, **Xb**, and **Xr** are defined for the frontal, rear, and roof recirculation of the building, respectively.

Table 3 presents a comparison of the length of different recirculation zones in comparison with available experimental data.


**Table 3.** Recirculation length of the recirculation zones.

**Figure 10.** Velocity streamlines of the mean flow field on the symmetry plane of the computational domain for Scenario A.

The recirculation zone **Xf** in front of the upstream surface of the building for Scenario A is 41% higher than the corresponding for Scenario B and 8.14% lower than Scenario C. In addition, the recirculation zone on the roof **Xr** for Scenario A is 20.81% lower than the recirculation zone of Scenario B and 27.44% lower than for Scenario C. The recirculation zone in the wake region **Xb** of the building for Scenario A is 13% lower than the corresponding recirculation zone of Scenario B and 23.28% lower than the recirculation zone for Scenario C.

Figure 11 shows the normalized vertical temperature profile for Scenario C at the position with coordinates of X = 7 H and Y = 5.5 H, on the symmetry plane and two heights behind the building.

The present numerical results are in good agreement with the experimental data of Uehara and Murakami [41], who studied the stability of the atmospheric flow inside an urban street canyon placed normal to the wind direction. They found that inside the canyon exists a stable thermal stratification as shown above, which weakens the cavity eddy.

Figure 12 shows the temperature distribution on the vertical and horizontal surfaces of the building for Scenario C, along the intersection lines with the transverse symmetry plane (path No. 3) and two other planes parallel to it at distances 0.25 H (path No. 2) and 0.45 H (path No. 1), respectively.

**Figure 11.** Normalized temperature profiles at the wake region of the flow for Scenario C, at X=7H and Y = 5.5 H, where *T* is the mean temperature, *Ta* is the ambient temperature, and *Tf* is the temperature at the floor of the computational domain [41].

The surface temperatures of the cube are the result of an energy balance between the incoming thermal radiation that is absorbed and reflected or emitted, and heat convection by the airflow. The amount of radiation that reaches the surfaces of the cube and the visibility of the solar path are determined by Athens's longitude and latitude. Albedo and emissivity of the surfaces control the shortwave and longwave radiation components reflected and emitted, respectively. The albedo and emissivity coefficients are 0.2 and 0.95, respectively, typical values for concrete material. The thermal heat capacity and thermal conductivity of the cube control its ability to store and conduct heat, respectively.

**Figure 12.** Temperature distribution on the vertical and horizontal surfaces of the building for the Scenario C along three different streamwise paths.

The temperature on the cube surfaces depends primarily on the convective thermal energy loss. In regions, where turbulence exists high heat exchange takes place. At the lower part of the frontal surface of the cube, where the horseshoe-shaped vortex appears, low surface temperatures appear. On the contrary, on the upper part of the frontal surface of the cube higher surface temperatures are shown. Similarly, on the wake zone of the cube where the arc vortex is present, low surface temperatures are shown because of the high heat exchange rate. On the first part of the roof surface where the flow separates the heat exchange is small, so the surface temperature is high, and on the second part it is reduced where the flow reattaches.

Figure 13 present the vertical mean pollutant concentration profiles computed for the time period from 800–1400 (s), at which all variables are statistically stationary on the symmetry plane for both Scenarios A and C at the Cartesian coordinate positions of X: 7 H, for a, X: 7.2 H for b, and 7.4 H, for c. The height Z extends from 0 to 2 H.

In Figure 13, the highest mean concentration of the pollutant is observed for Scenario A, where a large amount of the pollutant is trapped inside the wake region. In contrast, in Scenario C where because of the high buoyancy forces the maximum pollutant concentration is small in the above region as the pollutant is shifted almost one cube height higher away from the floor, thus creating a recirculation region of larger length and height than that corresponding for Scenario A.

As the pollutant in Scenario A is trapped inside the wake region of the cubical building, it becomes very harmful to the pedestrian's health. In Scenario C, this problem is exempted to a large extent, as the high levels of pollutant concentrations are at much higher heights than those of the pedestrians. This behavior is the result primarily of the buoyant forces acting on the roof and on the leeward surfaces of the building.

Figure 14 present the distribution of the flux concentration on the symmetry plane because of the x-velocity component for Scenarios A and C at the x-positions: (a) (X: 7 H, Y: 5.5 H), (b) (X: 7.2 H, Y: 5.5 H), and (c) (X: 7.4 H, Y: 5.5 H) for heights Z that extend from 0 to 2 H.

**Figure 13.** *Cont*.

**Figure 13.** Vertical concentration profiles with the pollutant source at X = 7 H for Scenarios A and C at positions (**a**) (X: 7 H, Y: 5.5 H), (**b**) (X: 7.2 H, Y: 5.5 H), and (**c**) (X: 7.4 H, Y: 5.5 H).

In the wake region, up to the height of the cube, the flux concentration is similar for both cases A and C. At heights higher than the building, the concentration flux for Scenario A is larger of that of C because of the important influence of higher free stream air velocity. Additionally, the wake zone height in Scenario C exceeds the height of the building and the low values of air velocity produce low concentration fluxes.

Figure 15a–c present the distribution of the flux concentration on the symmetry plane because of the z-velocity component for Scenarios A and C at the x-positions: (a) (X: 7 H, Y: 5.5 H), (b) (X: 7.2 H, Y: 5.5 H), and (c) (X: 7.4 H, Y: 5.5 H) for heights Z that extend from 0 to 2 H.

**Figure 14.** *Cont*.

**Figure 14.** Concentration flux profiles that are due to the x−velocity component on the symmetry plane for Scenarios A and C at positions (**a**) (X: 7 H, Y: 5.5 H), (**b**) (X: 7.2 H, Y: 5.5 H), and (**c**) (X: 7.4 H, Y: 5.5 H) for heights Z that extend from 0 to 2 H.

**Figure 15.** *Cont*.

**Figure 15.** Distribution of the non-dimensional flux concentration at the z−component of the velocity on the symmetry plane for Scenarios A and C at positions (**a**) (X: 7 H, Y: 5.5 H), (**b**) (X: 7.2 H, Y: 5.5 H), and (**c**) (X: 7.4 H, Y: 5.5 H) for height Z that extends from 0 to 2 H.

For Scenario A, where the buoyancy forces do not affect the pollutant distribution, the concentration flux is toward the ground floor (Figure 15a,b) within the recirculation zone near the point of reattachment, in accord with the flow field of Figure 10, and away from the floor outside the recirculation zone. On the contrary, for Scenario C, the buoyancy forces lift the pollutants near positionX=7H (left end of source) where at positions X = 7.2 H and X = 7.2 H (center and right end of source, respectively) the buoyancy forces are counterbalanced by momentum flux forces. It is expected that by placing the pollutant source closer to the leeward face of the building, the lifting forces that act on the concentration would be higher, resulting in higher concentration fluxes away from the ground. Thus, a higher fraction of the pollutants will escape the recirculation zone.

### **7. Conclusions**

In the present work, three different scenarios are studied. The emphasis is on the flow and pollutant dispersion around an isolated cubical building in an open area and with a pollutant source placed on the ground and at a distance of one building height in the symmetry plane behind its rear surface. Scenario A with adiabatic walls of the building and without any buoyancy forces, Scenario B with the same constant temperature of 50 ◦C on all the surfaces of the building, and Scenario C with convective and radiative conditions imposed by the local microclimate.

Each scenario is simulated using the large-eddy simulation (LES) approach. The numerical results that describe the nature and the behavior of the flow around an isolated building are in good agreement with the corresponding experimental and numerical data that were used for their validation.

It is observed that at higher flow field temperatures, the length and the corresponding height of the recirculation are increased. Thus, the recirculation region for Scenario B has higher length than that for Scenario A and lower than that for Scenario C, for which exist the highest temperatures from the three scenarios studied.

The pollutant in Scenario A is trapped inside the wake region of the cubical building with its highest concentration levels at a height close to what a pedestrian would breathe, thus being very harmful to his health. In Scenario C, the higher temperatures together with the associated buoyancy forces that are developed raise most of the pollutants to a height greater than that of the building, allowing them to escape the recirculation region of the building and be carried away by the wind, without causing any serious harm to the health of a pedestrian.

It is expected that by placing the pollutant source closer to the leeward face of the building, the lifting forces that act on the concentration would be higher, resulting in higher concentration fluxes away from the ground. Thus, a higher fraction of the pollutants will escape the recirculation zone.

**Author Contributions:** Conceptualization, C.L.P. and K.V.; methodology, C.L.P., K.V., I.E.S. and I.C.L.; software, C.L.P. and A.V.P.; validation, C.L.P. and K.V.; formal analysis; C.L.P.; resources, I.E.S.; data curation, C.L.P.; writing-original draft preparation, C.L.P. and K.V.; writing-review and editing, I.C.L. and I.E.S.; visualization, C.L.P. and A.V.P.; supervision, K.V.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used in this study are available on request from the corresponding author. The data can be easily reproduced from the theoretical analysis described in the study.

**Acknowledgments:** This work was supported by computational time granted from the National Infrastructures for Research and Technology S.A. (GRNET S.A.) in the National HPC facility—ARIS under project ID pr011045- UrbanFirePlan2. The authors acknowledge and express particular appreciation to the PALM group at the Institute of Meteorology and Climatology (IMUK) of Leibniz Universität Hannover, Germany that has developed the Parallelized Large Eddy Simulation Model (PALM).

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

### **References**


13. Pisello, A.L.; Pignatta, G.; Castaldo, V.L.; Cotana, F. The Impact of Local Microclimate Boundary Conditions on Building Energy Performance. *Sustainability* **2015**, *7*, 9207–9230. [CrossRef]


### *Article* **Thermal Performance Enhancement Using Absorber Tube with Inner Helical Axial Fins in a Parabolic Trough Solar Collector**

**Mohammad Zaboli 1, Seyed Soheil Mousavi Ajarostaghi 2, Seyfolah Saedodin <sup>1</sup> and Mohsen Saffari Pour 3,\***


**Abstract:** In the present work, a parabolic trough solar (PTC) collector with inner helical axial fins as swirl generator or turbulator is considered and analyzed numerically. The three-dimensional numerical simulations have been done by finite volume method (FVM) using a commercial CFD code, ANSYS FLUENT 18.2. The spatial discretization of mass, momentum, energy equations, and turbulence kinetic energy has been obtained by a second-order upwind scheme. To compute gradients, Green-Gauss cell-based method has been employed. This work consists of two sections where, first, four various geometries are appraised, and in the following, the selected schematic of the collector from the previous part is selected, and four various pitches of inner helical fins including 250, 500, 750 and 1000 mm are studied. All the numerical results are obtained by utilizing the FVM. Results show that the thermal performance improvement by 23.1% could be achieved by using one of the proposed innovative parabolic trough solar collectors compare to the simple one. Additionally, the minimum and maximum thermal performance improvement (compare to the case without fins) belong to the case with P = 250 mm by 14.1% and, to the case with P = 1000 mm by 21.53%, respectively.

**Keywords:** solar energy; parabolic trough collector; thermal performance; turbulator; swirl generator; numerical simulation

### **1. Introduction**

Heat exchangers, which transfer thermal energy through direct and indirect contact between fluids, are considered an indispensable part of several industries, from pharmaceutical to petrochemicals. Indirect contact heat exchangers are extensively used in solar systems. Due to the increasing importance of solar energy, nowadays, improving the performance of solar systems is one of the most important challenges for human beings and researchers. One of the most efficient types of solar collectors is a parabolic trough solar collector (PTSC), which is employed in both domestic and power plant applications. Recently, many efforts have been made by scientists to increase the efficiency of this type of collector. Actually, the PTSC component is a heat exchanger in which heat transfer fluid (HTF) flows in the receiver tube and absorbs the radiated solar energy.

In order to improve the thermal performance of this type of heat exchanger, various methods have been proposed and studied by researchers. Berger et al. [1] categorized heat exchangers' augmentation procedures into passive and active subdivisions. One of two categories of ameliorating heat exchange methods is passive procedures [2]. It means that there is no requirement for any type of additional force. This method includes techniques such as nano particle [3], helical tubes [4], treated area [5], vortex generator [6], displaced increase devices [7], extended surfaces [8], and jagged surfaces [9].

**Citation:** Zaboli, M.; Mousavi Ajarostaghi, S.S.; Saedodin, S.; Saffari Pour, M. Thermal Performance Enhancement Using Absorber Tube with Inner Helical Axial Fins in a Parabolic Trough Solar Collector. *Appl. Sci.* **2021**, *11*, 7423. https:// doi.org/10.3390/app11167423

Academic Editor: Rüdiger Schwarze

Received: 20 May 2021 Accepted: 2 August 2021 Published: 12 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/).

Tang et al. [10] checked the heat transfer in the heat exchanger. According to their results, an increase of length enhances the performance of the vortex-generator fin. Azari et al. [11] checked the heat exchange in a heterogeneous circular microchannel, the results show the maximum heat exchange rates are got when the boundary conditions are symmetrical. Darzi et al. [12] have performed an experimental probing in a corrugated tube with nanofluid, they deduced that nanofluid with corrugated tubing could increase heat exchange by 330% compared to net liquid. Du et al. [13] numerically have worked on tube layouts, they demonstrated that the comprehensive efficiency of the optimized heat transfer device is ameliorated by about 50%. Bahiraei et al. [14] have checked on a triple-tube that their outcomes show that, based on heat exchange efficiency, an increase in height and reduced pitch are suggested. Abolarin et al. [15] have checked the impact of various geometries in twisted tapes, and the results pointed out that connection angle increases heat transfer. They found that when the wavelength is smaller, the temperatures near the walls increased. Zhang et al. [16] checked the thermal performance in corrugated pipes. According to the outcomes, the corrugated pipe ameliorated the heat transfer rate. Therefore, the other passive methods containing turbulator, etc., are usually better in the turbulent flow [17]. Several of the new works related to the utilization of inserts and the other passive methods to improve the thermal performance of the heat exchangers are displayed and listed in Table 1.


**Table 1.** Overview of articles to improve the heat transfer.


**Table 1.** *Cont.*

Solar energy is considered as a heat flux for solar collectors. The survey of a solar heat transfer device is trusted for increasing the thermal performance of the available system to obtain an outcome [41]. Moghaddaszadeh et al. [42] used two passive techniques to ameliorate the efficiency of heat exchangers in solar collectors. The outcomes illustrate that the use of nanofluid went up the Nusselt number up to 4%. Ghasemi and Ranjbar [43] in the numerical verification on solar parabolic collector demonstrated that adding nanoparticles increase thermal performance and using Al2O3 or CuO as nanofluid at 3% volume fraction enhances the heat exchange up to 28 and 35 percent, respectively. Reddy et al. [44] investigated the thermal analysis of parabolic trough collectors for various geometrical parameters such as porosity and fin aspect ratio. Their results pointed out that the porosity will increase the heat exchange by about 17%.

According to the presented literature review of previous works (Table 1) related to heat transfer enhancement of the PTSC by various passive methods, it can be concluded that among different kinds of utilized methods, employing helical axial fins has not been considered and analyzed previously. In the present work, one parabolic trough solar collector with inner helical axial fins as swirl generator or turbulator is considered and analyzed. This work consists of two sections. In the first section, four various schematics of the collector are appraised. In the second section, the selected geometry of the collector from the previous part (first section) is selected, and four various pitches of inner helical fins including 250, 500, 750 and 1000 mm are studied.

#### **2. Materials, Methods and Boundary Conditions**

Parabolic trough collector systems, as demonstrated in Figure 1a, in the studied geometry, a parabolic concentrator is used, which is due to the reduction of heat loss and increase of the produced temperature. Figure 1 shows the schematic of the checked geometry, and also Table 2 displays the geometrical constants and parameters. According to the check accomplished, for turbulent flow the entry region length is five times more than the diameter of the inner tube [35]; however, some researchers are of the opinion that it has a slight impact [45,46]. The three-dimensional numerical simulations have been done by finite volume method using commercial CFD code, ANSYS FLUENT 18.2. The geometrical parameters of the considered heat exchanger are illustrated in Figure 1c. The <sup>30</sup> ◦C inflow enters the tube at different velocities, namely 0.2, 0.24, 0.28 and 0.32 m·s−1. Velocity–inlet and pressure–outlet boundary conditions were considered at the inlet and outlet ports, respectively. A heat flux of 60,000 W·m−<sup>2</sup> (is calculated by dividing the power consumption (6000 W) to the side surface area (0.10 m2) is transferred and applied to the outer part of the geometry made of an internally threaded steel layer. Moreover, the thermos–physical properties of tube material and fluid are noted in Table 3. The following situations are assumed in simulations:



**Table 2.** Physical parameters of the analyzed collector.

**Table 3.** The thermos–physical properties of materials [46].


The present work includes two sections. In the first part, three geometries with various schematics of the inner helical axial fins are considered and investigated, and the obtained results are compared with the simple collector without any insert (fins). The schematics of the evaluated geometries in the first section of present work are illustrated in Figure 2. In the second section, the selected model according to the first section obtained results is considered with various pitch of the inner helical axial fins (P) including 250, 500, 750 and 1000 mm which the evaluated cases are depicted in Figure 3.

**Figure 1.** (**a**) A schematic of the PTC structure, (**b**) the schematic of the considered pipe and the inner helical axial fins, and (**c**) the geometrical parameters.

**Figure 2.** The schematics of the evaluated geometries in the first section of present work.

**Figure 3.** The schematics of the evaluated geometries in the second section of present work. (**a**) P = 250 mm, (**b**) P = 500 mm, (**c**) P = 750 mm and (**d**) P = 1000 mm.

### **3. Governing Equations and Dimensionless Parameters**

The actual system is governed by the continuity, momentum, and energy equations that are defined, respectively, as [47,48]:

$$\frac{\partial \rho}{\partial t} + \nabla \nabla \cdot \left(\rho \stackrel{\rightarrow}{\boldsymbol{v}}\right) = \mathcal{S}\_{\text{ll}} \tag{1}$$

$$\frac{\partial \left(\rho \stackrel{\rightarrow}{\upsilon}\right)}{\partial t} + \nabla \cdot \left(\rho \stackrel{\rightarrow}{\upsilon} \stackrel{\rightarrow}{\upsilon}\right) = -\nabla p + \nabla \left(\stackrel{\rightarrow}{\tau}\right) + \rho \stackrel{\rightarrow}{\mathcal{g}} + \stackrel{\rightarrow}{F} \tag{2}$$

$$\frac{\partial(\rho E)}{\partial t} + \nabla \cdot \left(\vec{v}\left(\rho E + p\right)\right) = \nabla \left(k\_{eff}\nabla T - \sum\_{j} h\_{j}\vec{\tilde{J}}\_{j} + \left(\overline{\tilde{\pi}}\_{eff}, \overrightarrow{\tilde{\upsilon}}\right)\right) + S\_{h} \tag{3}$$

The thermal performance coefficient (η) can be exerted to quantify the yield of a heat exchanger that contains the Nusselt number and the friction factor (*f*), and is defined as [49–58]:

$$N\mu = \frac{h\_m d\_h}{k} \tag{4}$$

$$f = \frac{2d\_h \Delta P}{\rho u^2 L} \tag{5}$$

$$\eta = \left(\frac{Nu}{Nu\_0}\right) \left(\frac{f\_0}{f}\right)^{\frac{1}{3}} \tag{6}$$

The three-dimensional numerical simulations have been done by finite volume method using commercial CFD code, ANSYS FLUENT 18.2. The spatial discretization of the mass, momentum, turbulence kinetic energy, turbulence dissipation rate and energy equations has been achieved by a second-order upwind scheme. The velocity-pressure coupling has been overcome by the SIMPLE algorithm. To calculate the gradients, Green-Gauss cell-based method has been utilized. The convergence criteria were set to 10−<sup>6</sup> for the residuals of all equations except energy equation (10<sup>−</sup>8).

### **4. Results**

### *4.1. Grid Independency Study*

The generated grids for the proposed collector geometry is displayed in Figure 4. Four variant grids were checked to certify the accuracy of the numerical outcomes. The analysis outcomes are offered in Table 4. It is vivid that the two middle mesh sizes do not vary substantially and in the interest of minimizing the computational time whilst attaining trust outcomes, demonstrate that a third option is a good selection.

**Figure 4.** Considered grid for proposed geometry as PTC.



### *4.2. Validation Study*

The Dittus-Boelter correlation [59] provides the Nusselt number as a function of the Re and Pr according to Equation (7):

$$\text{Nu} = 0.023 \,\text{Re}^{0.8} \,\text{Pr}^{0.4} \tag{7}$$

Additionally, the Seban-Shimazaki correlation [59] (Equation (8)) is employed to compare with the simulation outcomes. A high correlation is suggested for PeD ≥ 100. The outcomes of the validation study are demonstrated in Figure 5. Therefore, the available work presents good accuracy. Additionally, to validate the present numerical model with the results of the performance of a real PTC system, the experimental results of Zou et al. [60] are used, and the comparison is presented in Figure 6. It can be concluded that present numerical method has good capacity for modeling the fluid and thermal characteristics in a PTC system. The uncertainty in measurement of the experimental parameters of the Zou et al. work [60] was as follows. The uncertainty in HTF velocity measurement was ±0.05m/s for HTF velocity up to 0.28 m/s and ±0.2 m/s for HTF velocity between 0.28 and 0.32 m/s. The uncertainty of HTF temperature was ±0.3◦C in this study.

$$\text{Nu} = \text{5} + 0.025 \,\text{Pe}\_{\text{D}}^{0.8} \tag{8}$$

**Figure 5.** Validation betwixt the available simulation and (**a**) Dittus–Boelter correlation [41]; (**b**) Seban–Shimazaki correlation [59].

**Figure 6.** Results of the validation analysis between the present numerical work and experimental work of Zou et al. [60].

### *4.3. The Impact of the Schematic of the Proposed Collector's Fins*

In this section, the influence of the collector shape on thermal efficiency is investigated numerically. Three various models are considered, and the obtained results have been compared with the simple shaped collector. The considered models here are shown in Figure 2. The variation of heat transfer coefficient versus different inlet velocity for various models is illustrated in Figure 7. Hence, it can be concluded that firstly, all the proposed models show higher heat transfer coefficient more than the simple collector (Model 1). Secondly, it is depicted that among the studied models, the maximum heat transfer coefficient belongs to Model 2 in all studied inlet velocities. Additionally, the minimum heat transfer coefficient belongs to Model 4 in all considered inlet velocities.

**Figure 7.** Heat transfer coefficient variation with the inlet velocity for various models at P = 500 mm.

To realize better the effect of geometry on the thermal performance of the proposed system, the streamline with the contour of temperature for all four studied models is shown in Figure 8. Hence, it can result that the schematic of the geometry and also the attendance of the inner axial helical fins have a significant impression on the rotation flow production and finally more heat transfer rate. Furthermore, contours show that Model 4 could enhance the heat transfer rate more than the other models. Additionally, all models depict better thermal performance in comparison with the collector equipped with simple channel (Model 1) because of the presence of swirl flows generated by the helical fins.

**Figure 8.** The streamline with contour of temperature for various models at P = 500 mm.

The variation of pressure drop and friction factor (f) versus various inlet velocities for various models are illustrated in Figure 9a,b, respectively. Figure 9a depicts that all models have more pressure drop compare to Model 1 (collector with simple channel, without helical fins) obviously because of presence of swirl flows. The maximum pressure drop belongs to Model 2 and the Model 4 shows lowest pressure drop among these four models. Additionally, as the inlet velocities rises, the pressure drop increases due to the effect of forced convection in the channel. This trend is the same for all models. Moreover, by increasing the inlet velocities, the differences between the studied models rises. This point could be important considering the pressure drop factor, as the differences between Model 4 and Model 1 (as base model) are very low.

**Figure 9.** (**a**) Pressure drop and (**b**) friction factor variation with the inlet velocity for various models at P = 500 mm.

Figure 9b shows that the trend of variation of fraction factor among the models is the same with pressure drop (Figure 9a). However, it should be focused that the variation of friction factor with inlet velocities is not the same with pressure drop, and they are exactly the opposite. By increasing the inlet velocity, the friction factor declines. According to Equation (5), the friction factor is directly related to the pressure drop; however, the squared velocity is placed at the denominator of the equation.

Contours of temperature at outlet for various models and P = 500 mm are presented in Figure 10. The impact of changing the schematic of the collector on temperature distribution in the outlet port and also the heat transfer rate is presented clearly in this figure. Accordingly, all models show better temperature distribution than the Model 1 (simple one, without fins).

**Figure 10.** Contours of temperature at outlet and a slice (X = 0) for various models at P = 500 mm and VInlet = 0.1256 m/s.

The best parameter to comprehensively study the impact of efficient parameters on the performance of the proposed collector is thermal efficiency (η), which is calculated by Equation (6). In this Equation, the index 0 refers to the base collector (without fins). The variation of this parameter versus various inlet velocity for variant models are presented in Figure 11. Hence, it can be firstly concluded that all the models illustrate higher thermal performance than Model 1 (simple collector without fins), which means that the values of thermal performance are higher than unity. Furthermore, the highest and lowest thermal performance improvements (compared to Model 1) belong to the Model 3 by 23.1% (at Vinlet = 0.2512 m/s) and Model 3 by 20.7% (at Vinlet = 0.1884 m/s), respectively. Contours of velocity magnitude at surface X = 0 for various models at P = 500 mm and VInlet = 0.1256 m/s are illustrated in Figure 12. It can be realized that among the considered models in this section, Model 2 depicts higher velocity magnitude and consequently more swirl flows in the proposed collector, which leads to a higher heat transfer rate.

**Figure 11.** The thermal performance (η) versus inlet velocity for various models.

**Figure 12.** Contours of velocity magnitude at surface X = 0 for various models at P = 500 mm and VInlet = 0.1256 m/s.

### *4.4. The Impact of the Inner Helical Fins Pitch (P)*

In the second section, according to the obtained numerical results from the last section (Section 4.3), Model 2 (see Figure 2) is considered as the selected model here. Results from the last section showed that Model 2 has the maximum heat transfer coefficient with the highest pressure drop. Therefore, to study the impact of the inner helical fins pitch, this model is used. In this section, the impact of the inner helical fins pitch on the thermal efficiency of the system is evaluated numerically. Four various pitches of the inner helical fins including 250, 500, 750 and 1000 mm are considered, and the obtained outcomes are compared with the simple collector (without fins). The considered models here are shown in Figure 3. The variation of heat transfer coefficient versus different inlet velocity for various models is illustrated in Figure 13.

**Figure 13.** Heat transfer coefficient variation with the inlet velocity for various pitches of inner helical axial fins.

According to Figure 13, it can be concluded that firstly, all the proposed models show higher heat transfer coefficient more than the simple collector (without fins). Secondly, it is depicted that among the studied models, the unmost heat transfer coefficient belongs to case with P = 250 mm in all studied inlet velocities. Additionally, the least heat transfer coefficient belongs to the case with P = 1000 mm in all considered inlet velocities. To realize better the impact of the pitch of inner helical fins on the thermal performance of the proposed collector, the contour of temperature at outlet and a slice in the middle of the computational domain (X = 0) for all cases are shown in Figure 14. It can be obviously seen that higher temperature in the proposed collector is achieved at lower pitch of the inner helical fins.

**Figure 14.** Contours of temperature at outlet and a slice (X = 0) for various models at VInlet = 0.1256 m/s.

The variations of pressure drop and friction factor (f) versus various inlet velocities for various models are illustrated in Figure 15a,b, respectively. Figure 15a depicts that all models have higher pressure drop compared to the case without helical fins obviously because of the presence of swirl flows. The greatest pressure drop belongs to the case with P = 250 mm, and the case with P = 1000 mm shows the lowest pressure drop among the investigated models. Additionally, as the inlet velocities rise, the pressure drop increases due to the effect of forced convection in the channel. This trend is the same for all models. Moreover, by increasing the inlet velocities, the differences between the studied models

rise. This point could be important considering the pressure drop factor, as the differences between the cases with P = 500, 750 and 1000 mm are very low.

**Figure 15.** (**a**) Pressure drop and (**b**) friction factor (f) variation with the inlet velocity for various models at P = 500 mm.

Figure 15b shows that the trend of variation of fraction factor among the models is the same with pressure drop (Figure 15a). On the other hand, as same as Figure 15a, the highest and lowest friction factor belong to the cases with P = 250 mm and P = 1000 mm, respectively. However, it should be focused that the variation of friction factor with inlet velocities is not the same with pressure drop and they are exactly the opposite. By increasing the inlet velocity, the friction factor declines. According to Equation (5), the

friction factor is directly related to pressure drop; however, the squared velocity is placed at the denominator of the equation.

The streamline with the contour of temperature for various models are presented in Figure 16. The impact of changing the schematic of the collector on temperature distribution in the outlet port and also the heat transfer rate is presented clearly in this figure. Accordingly, all models show better temperature distribution than the case without fins.

**Figure 16.** *Cont*.

**Figure 16.** The streamline with contours of temperature for various models.

The best parameter to study comprehensively the impact of efficient parameters on the performance of the proposed collector is thermal efficiency (η), which is calculated by Equation (6). The variation of this parameter is presented in Figure 17. Therefore, it can be concluded firstly that all the models illustrate higher thermal performance than the case without fins, which means that the values of thermal performance are higher than unity. Furthermore, the greatest and smallest thermal performance improvement (compare to the case without fins) belong to the case with P = 1000 mm by 21.53% and the case with P = 250 mm by 14.1% for the latest velocity, respectively.

**Figure 17.** The thermal performance (η) versus various inlet velocities for different models.

### **5. Conclusions and Future Scope**

In the present work, a PTC collector with inner helical axial fins as a swirl generator or turbulator is considered and analyzed. This work consists of two parts. In the first part, four various geometries of the collector are appraised. In the second part, the selected schematic of the collector from the first sector is selected and four various pitches of inner helical fins including 250, 500, 750 and 1000 mm are studied. Results presented that the thermal performance improvement by 23.1% could be achieved by using one of the proposed innovative parabolic trough solar collector compare to the simple one. Additionally, the utmost and least thermal performance improvement (compare to the case without fins) belong to the case with P = 1000 mm by 21.53% (at Vinlet = 0.314 m/s) and the case with P = 250 mm by 14.1% (at Vinlet = 0.314 m/s), respectively.

Although present work has presented new insight into the flow and thermal characteristics of a parabolic trough solar collector combining innovative inner helical axial fins as a swirl generator, additional research is required to determine the capability of this design in more realistic engineering situations by the development of an experimental prototype.

**Author Contributions:** Investigation, Software, Validation, Writing—original draft, M.Z.; Formal analysis, Methodology, Validation, Writing—review & editing, S.S.M.A.; Supervision, Writing review & editing, S.S.; Funding acquisition, Writing—review & editing, M.S.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data available on request due to restrictions eg privacy or ethical.

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

**Expects Data:** Data available on request from the authors: The data that support the findings of this study are available from the corresponding author upon reasonable request.

### **Nomenclature**



### **References**


### *Article* **Experimental and CFD Investigation into Using Inverted U-Tube for Gas Entrainment**

### **Khaled Yousef 1, Ahmed Hegazy <sup>1</sup> and Abraham Engeda 2,\***


Received: 3 November 2020; Accepted: 14 December 2020; Published: 18 December 2020

**Abstract:** An experimental and numerical study is presented in the current work for gas entrainment using an inverted vertical U-tube. Water flows vertically up in an inverted U-tube which creates a low-pressure region in the tube upper portion. This low-pressure region can be used to extract gases by connecting it to a branch pipe. The extracted gases considered in this work are a mixture of air and water vapor. The water vapor from the side branch pipe is mixed with the flowing water under the siphon effect. This results in a progressive water vapor condensation as the mixture proceeds towards the exit due to an increase in vapor partial pressure. The air is drawn by inertia to be released out at the tube lower exit of the inverted U-pipe. The current study deals with these complicated flow behaviors due to the mixing undergoing condensation. A test rig is designed for experimentally studying the behavior of water flow in an inverted U-tube where the air is mixed with the flowing water at the top region of this tube. The CFD computations are accomplished for a side gas mixture with volume fractions up to 0.7 with water vapor mass fractions in this mixture to be 0.1–0.5. The tested water mass flow rates in the main tube are 2, 4, 6, 8 kg/s to account for all possible flow mass ratios. The CFD computations are validated with water and air two phase flow with the measurements of both the experiments of the current research and the literature. The present results reveal that slightly raising the water mass flow rate at a constant side mixture mass ratio produces a reduced generated pressure in the upper tube part. This is attributed to extra water vapor condensation taking place rapidly by increasing the water flow rate in the tube upper part. Furthermore, the turbulence quantities begin to break down at a side mixture volume fraction of 0.55 with water and air mass flow rates of 2 kg/s and 0.002 kg/s, respectively. On the other side, raising the air mass flow rate at the higher values of water vapor and water mass flow rates breaks the generated vacuum pressure and turbulence due to entrainment. Moreover, this proposed framework can produce a lower static pressure, reaching 55.1 kPa, which makes it attractive for gas extraction. This new technique presents innovative usage with less consumable energy for extracting gases in engineering equipment.

**Keywords:** inverted U-tube; two-phase flow; air/water vapor entrainment; steam condenser

### **1. Introduction**

Few engineering and medical applications need evacuated cavities during certain manufacturing and testing processes. These evacuated cavities may contain a variety of gases and liquids. So, the continuous removal of these undesired fluids and maintenance of the required evacuated cavities at a specific pressure becomes mandatory for these applications. The most used devices for venting these gases are categorized as a steam trap, liquid ring, vacuum pumps, manual valve, automatic valve and steam ejector. All these devices consume a considerable amount of energy to accomplish their duties. The use of an inverted U-tube is an innovative method for entraining such gases, which may lead to the consumption of less energy [1]. When water flows up in an inverted U-tube, it produces a low-pressure region at its upper part. This produces a low-pressure region, which can be exploited to draw gases from the required cavity by connecting the top of the inverted U-tube to the cavity through a side pipe. In the present paper, two-phase flow (water and gas) and three components (water liquid, air and water vapor) are considered in the proposed CFD computations. One of the most important pieces of equipment working under vacuum is the steam condenser of the steam power plant. The main function of the steam condenser, besides steam condensation, is to preserve the generated operational vacuum by evacuating the gases in the cycle. Currently, steam power plants share the largest percentage of electricity generation worldwide [2,3].

As commonly known in a condenser, the volume fractions of air, vapor, and water change rapidly, and this complicated behavior has a significant impact on plant efficiency. Water vapor can be taken out of the condenser cavity through condensation and it is pumped back to the water feeding system of the plant boiler. The non-condensable gas is mainly air that has leaked into the plant turbine and condenser, which are working below the atmospheric pressure or it comprises air and non-condensable gases released from the plant working fluid. Furthermore, these non-condensable gases may also be formed by water decomposition. For efficient plant operation, these gases are extracted out of the condenser cavity to avoid increasing the condenser operating pressure. Raising the condenser operating pressure diminishes the steam turbine output power. Additionally, these gases severely reduce the heat transfer rates in the condenser by the covering wrap of its tubes. Furthermore, the oxygen content raises the tube corrosiveness in the condenser, which lowers the condenser's operational life. So, these gases should be evacuated continuously for better operation and longer life of the steam condenser [4–6].

Few experiments consider the two-phase flow of steam condensation with the presence of gas in vertical tubes [7–9]. The main finding of these experiments was that raising both the gas mass fraction and the inlet properties reduced the heat transfer rate. Additionally, the local heat transfer was measured in the case of steam only and steam mixed with helium and air. The obtained temperature profiles for the considered cases provided a reasonable way of validating the proposed numerical model. The CFD computations of condensation in the existence of gases encountered some implementing difficulties but the resulted CFD flow profiles helped to investigate the physical mechanisms due to condensation. The effect of non-condensing gases on the steam condenser performance is presented theoretically by Strušnik et al. [10]. In their study, the steam ejector pump system is modeled by monitoring the gas extraction. Optimizing the steam ejector pump geometry was recommended for enhancing plant efficiency [10]. Strušnik et al. [11] made a comparative study for the power consumptions in the case of using the ejector and electric vacuum pump. They presented an economic analysis and operating costs guidelines for the appropriate selection of the condenser evacuating system in the steam power plant. It was found that the considered ejector is more consuming for the plant energy, as compared to the liquid ring vacuum pump. The condensation of steam with the existence of non-condensing gas was studied in [12,13]. They explored that the existence of the non-condensing gases severely reduced the heat transfer coefficient in the steam condenser. A hybrid system for gas extraction in the plant condenser was studied by Kapooria et al. [4]. In this hybrid system, a steam ejector integrated with a liquid ring vacuum pump was utilized. This proposed hybrid system was able to generate more vacuum in the steam condenser, while an observed increase in energy consumption was recorded.

Significant attention has been directed toward the complicated air and water flow patterns and their application in the steam condenser. In this regard, Yao et al. [14] studied numerically the two-phase flow of air/water in a small horizontal pipe. The obtained CFD results had an acceptable agreement with the considered experimental data and the CFD simulation accurately predicted the flow structures. The noticed small discrepancies between the experimental data and the CFD results may be attributed to the requirement for the full three-dimensional CFD simulations. In the same way, Juggurnath et al. [15] numerically presented the various two-phase flow profiles created in an inclined tube. They considered the effect of both surface tension and gravitational forces on the generated flow patterns. The air/water flow in a horizontal pipe was investigated numerically by Vásquez et al. [16]. Different multiphase flow models were used for obtaining the flow patterns and the proposed CFD model was validated. The validation showed a reasonable agreement with slight deviations. Seven experimental facilities for the air and water co-current flow in sloping pipes were presented by Pothof and Clemens [17]. A CFD geometrical and operational analysis for preventing the air accumulation in the tube was proposed. Panella [18] numerically and experimentally predicted the air and water flow in pipes. In this study, a framework for mixing air and water in the tube is fabricated. A circulating pump was utilized for controlling the water mass flow rate. This study revealed good validation in a wide range of data. Air and water two-phase flow in a tube with gas separator was presented by Afolabi and Lee [19]. This study ensured the ability of CFD models to accurately predict the behavior of the two-phase flow. Besides, the computational results had a good validation with the experiments.

The objective of the present work is to numerically capture the two-phase flow patterns in the inverted U-tube when mixing water with a side mixture of air and water vapor. Furthermore, the required flow conditions and operational parameters for using the inverted U-tube for the side entrainment of the air-water vapor mixture will be determined. Additionally, a reduced scale experimental test rig is fabricated to validate the present CFD results.

### **2. Experimental Setup**

An experimental setup is fabricated to disclose the siphon flow characteristics with side air entrainment in the inverted U-tube. The experimental test rig with half scaled dimensions, compared to the CFD model dimensions, is fabricated. The experimental set up comprises two tanks with a water elevation difference of ΔZ—see Figure 1. An inverted U-tube of 2.54 cm inner diameter, with its inlets immersed into the two water tanks, is utilized. The water levels in the tanks are changeable by lowering/lifting the right-side tank. The main filling port (2) in the left side tank and a water exit port in the right-side tank (12) is used to maintain the water levels between the two tanks constant during measurements. Bourdon pressure gauges (4) are utilized to measure the flow pressure around the system. The flowing water through the inverted U-tube is measured by using a calibrated tank and stopwatch while the side air is measured by using a variable area rotameter (7). The rotameter has a full-scale accuracy/repeatability of 6/2%. The rotameter has a variable flow area valve (8) to control the inlet air flow rate (9). The siphon flow generation is one of the tricky points before measurements. The pre-filling port (5) is used to fill the whole system with water while the air venting port (6) is used to release the air from the system cavity. The inverted U-tube inlets in the right- and left-side tanks are blocked during this system pre-filling by using manual valves.

Once the whole system is full of water and free of air, the inverted U-tube inlets in the right- and left-side tanks are unblocked while maintaining the water pre-filling port (5) opened to help initiate the siphon flow for a while. Reaching a steady siphon flow is the main task that should be performed first. Then, the desired water elevation difference between the left and right tanks is obtained and kept constant. The water flow rate is measured before and after side air enters the system at the same ΔZ. The air flow rate is changeable by using the variable area valve (8), which is also used to close the side pipe during the siphon flow initiation. The water level difference between the two tanks can be changed by adjusting the side air entrainment and the experiments are repeated for different flow conditions.

**Figure 1.** Experimental test rig and its components.

### **3. Computational Model**

The two-phase model is a set of equations and mathematical relations that are used to describe the phase's motion and interactions. The most commonly utilized two-phase models in the literature are homogeneous, drift flux and separated flow models [20–25]. The homogeneous model is less complicated and can be applied to many complex two-phase flows with good accuracy. The homogeneous model is represented by the mixture two-phase model in CFD computations, as will be discussed later in this paper. In homogenous models, the average properties for the flow variables are considered in the calculations and the relative motion between phases can be accurately calculated [20,26]. In addition, the homogenous model is easy to implement as its governing equations are close to the single-phase flow equations. In the present study, the mixture model for calculating the two-phase flow variables in the inverted U-tube is utilized.

Figure 2 shows the CFD model with its boundary conditions. In the proposed work, the siphon flow generated by water flowing up in an inverted U-tube creates a vacuum in the tube upper region. This vacuum region can be connected by a side pipe to extract a mixture of air and water vapor from the equipment cavity, the equipment cavity is not shown in this figure. The equipment cavity is considered the steam condenser in the present case. This results in mixing the water with the induced

side mixture of air and water vapor. This mixing will cause the water vapor condensation while all the flow mixture will continue flowing down due to flow inertia and gravity. At the lower outlet of the inverted U-tube, the air is released to the environment while water is re-circulated through the inverted U-tube. Furthermore, a variable speed circulating pump is used to transfer the water between the two tanks at the inlet and outlet of the inverted U-tube (the pump is not shown in the figure). Additionally, the hydrostatic head difference between the water surfaces in both tanks can be altered to control the water mass flow rate. The considered tube has low maintenance and operational cost, which makes it more favorable compared to the other energy consuming systems, i.e., water ring pump and steam ejector [27].

**Figure 2.** CFD model for the inverted U-tube with dimensions.

### *3.1. Governing Equations*

Steady-state three dimensional CFD computations were performed for the considered geometry shown in Figure 2. Differential equations governing the fluid flow are given as [28,29].

$$
\Delta \left( \rho \mathbf{u} \right) \tag{1}
$$

$$\Delta \left( \rho \mathbf{u} \,\mathbf{u} \right) = -\frac{\partial \mathbf{p}}{\partial \mathbf{x}} + \Delta \left( \mu \,\text{grad } \mathbf{u} \right) + \left[ -\frac{\partial \overline{\rho \mathbf{u} \mathbf{u}'^2}}{\partial \mathbf{x}} - \frac{\partial \overline{\rho \mathbf{u} \mathbf{u}' \mathbf{v}'}}{\partial \mathbf{y}} - \frac{\partial \overline{\rho \mathbf{u} \mathbf{u}' \mathbf{w}'}}{\partial \mathbf{z}} \right] \tag{2}$$

$$\Delta\left(\text{pu }\mathbf{v}\right) = -\frac{\partial\mathbf{p}}{\partial\mathbf{y}} + \Delta\left(\mu\,\text{grad }\mathbf{v}\right) + \left[ -\frac{\partial\overline{\rho\,\mathbf{u'}\mathbf{v'}}}{\partial\mathbf{x}} - \frac{\overline{\partial\rho\,\mathbf{v'}^2}}{\partial\mathbf{y}} - \frac{\partial\overline{\rho\,\mathbf{v'}\mathbf{w'}}}{\partial\mathbf{z}} \right] \tag{3}$$

$$\Delta \left( \rho \mathbf{u} \,\mathrm{w} \right) = -\frac{\partial \mathbf{p}}{\partial \mathbf{z}} + \Delta \left( \mu \,\mathrm{grad} \,\mathrm{w} \right) + \left[ -\frac{\partial \overline{\rho \mathbf{u'} \,\mathrm{w'}}}{\partial \mathbf{x}} - \frac{\partial \overline{\rho \mathbf{v'} \,\mathrm{w'}}}{\partial \mathbf{y}} - \frac{\partial \overline{\rho \mathbf{w'^2}}}{\partial \mathbf{z}} \right] \tag{4}$$

The energy equation can be found as follows to consider the heat transfer, as in the current study.

$$
\Delta \left( \rho \text{Iu} \right) = \Delta \left( \text{k grad T} \right) + \Phi + \left[ -\frac{\overline{\partial \rho} \, \overline{\mathbf{u}' \, \mathbf{l}'}}{\partial \mathbf{x}} - \frac{\overline{\partial \rho} \, \overline{\mathbf{v}' \, \mathbf{l}'}}{\partial \mathbf{y}} - \frac{\overline{\partial \rho} \, \overline{\mathbf{w}' \, \mathbf{l}'}}{\partial \mathbf{z}} \right] \tag{5}
$$

I is the total energy, <sup>I</sup> <sup>=</sup> <sup>h</sup> <sup>−</sup> **<sup>p</sup>** <sup>ρ</sup> <sup>+</sup> (**u**2+**v**<sup>2</sup>) <sup>2</sup> , while <sup>Φ</sup> is the function of dissipation. <sup>h</sup> is taken to be the sensible enthalpy where h = i Yihi and hi <sup>=</sup> <sup>T</sup> Tref CpidT. Tref is the reference temperature and it is equal to 298.2 K.

### *3.2. Realizable <sup>k</sup>* <sup>−</sup> <sup>E</sup> *Model*

In the present study, the realizable *<sup>k</sup>* <sup>−</sup> <sup>ε</sup> model in the mixture form is utilized for computing the turbulence fields. The choice of turbulence model is very critical in any CFD simulation. The realizable *<sup>k</sup>* <sup>−</sup> <sup>ε</sup> model is can provide outstanding performance in many complicated flows. Additionally, this turbulent modeling is characterized over the other turbulence models as it has a modified equation for calculating turbulent viscosity and a new transport equation for the dissipation rate ε [28,29].

$$\frac{\partial}{\partial \mathbf{x}\_{\mathrm{i}}} (\rho \mathbf{k} \mathbf{u}\_{\mathrm{i}}) = \frac{\partial}{\partial \mathbf{x}\_{\mathrm{i}}} \bigg[ \left( \mu + \frac{\mu\_{\mathrm{t}}}{\sigma\_{\mathrm{k}}} \right) \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{\mathrm{i}}} \bigg] + \mathbf{G}\_{\mathrm{b}} + \mathbf{G}\_{\mathrm{k}} - \rho \boldsymbol{\varepsilon} - \mathbf{Y}\_{\mathrm{M}} + \mathbf{S}\_{\mathrm{k}} \tag{6}$$

$$\frac{\partial}{\partial \mathbf{x}\_{\mathrm{i}}} (\rho \varepsilon \mathbf{u}\_{\mathrm{i}}) = \frac{\partial}{\partial \mathbf{x}\_{\mathrm{i}}} \Big[ \Big( \mu + \frac{\mu\_{\mathrm{t}}}{\sigma\_{\mathrm{c}}} \Big) \frac{\partial \varepsilon}{\partial \mathbf{x}\_{\mathrm{i}}} \Big] + \rho \mathbf{C}\_{1} \mathbf{S} \varepsilon - \rho \mathbf{C}\_{2} \frac{\varepsilon^{2}}{\mathbf{k} + \sqrt{\nu \varepsilon}} + \mathbf{C}\_{1\varepsilon} \frac{\varepsilon}{\mathbf{k}} \mathbf{C}\_{3\varepsilon} \mathbf{G}\_{\mathrm{b}} + \mathbf{S}\_{\varepsilon} \tag{7}$$

In the mixture turbulence model, the same turbulence field is given for all phases. This means that all phase properties are replaced with mixture properties in the equations. Gk is the turbulent kinetic energy generation, due to the mean velocity gradients, while Gb is the turbulence kinetic energy generation due to buoyancy. Sk and S<sup>ε</sup> are user-defined source terms. The constants in these equations are [28,29].

$$\mathbf{C}\_{1} = \max\left[0.43, \frac{\eta}{\eta + 5}\right], \eta = \mathbf{S}\frac{\varepsilon}{\mathbf{k}}, \mathbf{S}\sqrt{2\mathbf{S}\_{\text{ij}}\mathbf{S}\_{\text{ij}}} \tag{8}$$

In turbulence modeling, <sup>μ</sup><sup>t</sup> <sup>=</sup> <sup>ρ</sup> <sup>C</sup><sup>μ</sup> <sup>k</sup><sup>2</sup> <sup>ε</sup> , μeff = μ + μt.

The favorable features of the considered turbulence model, as discussed earlier, is Cμ, which does not have a fixed value and can be correlated as [28–30].

$$\mathbf{C}\_{\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\tiny\text{\text}}}}}}}}}}}}}}}\mathbf}}\mathbf} \text{\text{\text{\tiny\text{\text{\tiny\text{\tiny\text{\textdegree}}}}}\text{U} = \sqrt{\text{S}\_{\text{\text{\tiny\text{\textgreater\text{\tiny\text{\textgreater}}}}} + \hat{\Omega}\_{\text{\text{\textquad}}}} \hat{\Omega}\_{\text{\text{\textquad}} \text{\text{\textquotedblleft}}} \,} \,\,\hat{\Omega}\_{\text{\textquotedblright}} = \,\,\overline{\Omega}\_{\text{\textquotedblleft}} - \varepsilon\_{\text{\textquotedblleft}} \,\omega\_{\text{\textquotedblleft}} \,\omega\_{\text{\textquotedblright}}} \tag{9}$$

Here, <sup>Ω</sup>ij is the rotation tensor rate with angular velocity <sup>ω</sup><sup>k</sup> and A0 <sup>=</sup> 4.04, As <sup>=</sup> <sup>√</sup> 6 cosϕ, ϕ = <sup>1</sup> <sup>3</sup> cos−<sup>1</sup> √ 6*W* , *<sup>W</sup>* <sup>=</sup> SijSjkSki S <sup>3</sup> , S = SijSij, Sij = <sup>1</sup> 2 ∂Ui ∂xj + ∂Uj ∂xi . In the current study, C1<sup>ε</sup> = 1.44, C2 = 1.9, σ<sup>k</sup> = 1.0, and σ<sup>ε</sup> = 1.2.

### *3.3. Two-Phase Model*

The mixture two-phase model in Ansys Fluent is used for the present computations [31]. This model is the best representation of the homogenous two-phase flow for water, air and water vapor flow in the U-tube. In the mixture model, all phases share the same cell volume and can interchange energy and momentum. The cell volume fraction for the phase is α, and all phase volume fraction summation is equal to 1. In the following equations, α<sup>m</sup> is the volume ratio of the side mixture in a combination of the mixture (water vapor and air) and water flows. Xw and Xv are the water and water vapor mass flow rate ratios concerning the side mixture mass flow rate, respectively. Ya is the ratio of air mass flow rate to the total fluid mass flow rate, while Va,v is the volume for the side mixture [29,32,33].

$$\sum\_{\mathbf{j=1}}^{n} \alpha\_{\mathbf{j}} = 1$$

$$\alpha\_{\rm m} = \frac{\rm V\_{\rm a,V}}{(V\_{\rm w} + V\_{\rm a,v})}, \ \dot{m}\_{\rm a,v} = \dot{m}\_{\rm a} + \dot{m}\_{\rm v} \tag{11}$$

$$\mathbf{X\_{w}} = \frac{\dot{\mathbf{m\_{W}}}}{\dot{\mathbf{m\_{3,V}}}}, \mathbf{X\_{v}} = \frac{\dot{\mathbf{m\_{V}}}}{\dot{\mathbf{m\_{3,V}}}}, \mathbf{Y\_{a}} = \frac{\dot{\mathbf{m\_{3}}}}{\dot{\mathbf{m\_{W}}} + \dot{\mathbf{m\_{3}}}} \tag{12}$$

For the filled cells with water in two-phase flows, the volume fraction of water, αw, is 1 while α<sup>v</sup> and α<sup>a</sup> are taken to be 0 and the same way for the filled water vapor cells. If the cell combines the mixing of air, water, and water vapor, then these cell volume fractions are between 0 and 1. For calculating the secondary and primary phase volume fractions, FLUENT® solves separate continuity equations for secondary phases, while Equation (10) is utilized for computing the primary phase volume fraction. In the condensation process, the water is generated from water vapor due to phase change. Therefore, the water is taken to be the primary phase while water vapor and air are the secondary phases. The continuity equation for calculating the secondary phase volume fraction is given as follows [29].

$$\nabla \cdot \mathbf{x}\_{\text{scc}} \rho\_{\text{scc}} \stackrel{\rightarrow}{\mathbf{v}}\_{\text{m}} = \sum\_{j=1}^{n} \dot{\mathbf{m}}\_{pri-\text{scc}} - \dot{\mathbf{m}}\_{\text{scc}} \cdot \text{pri} \tag{13}$$

. <sup>m</sup>*pri*−sec, . msec−*pri* is a mass source term that represents the mass transfer rate from primary phase to secondary phase in the case of evaporation or mass transfer rate from secondary phase to primary phase in the case of condensation, respectively. The average velocity of the mixture is, <sup>→</sup> vm = *n <sup>j</sup>*=<sup>1</sup> αjρ<sup>j</sup> → vj ρ<sup>m</sup> and the mixture density is computed for all phase densities, <sup>ρ</sup><sup>m</sup> <sup>=</sup> *<sup>n</sup> j*=1 αjρj. All the mixture properties are correlated, likewise, with these mixture density calculations.

The source term added in the energy equation and the internal energy is found as

$$\mathbf{S}\_{\rm I} = \left(\dot{\mathbf{m}}\_{\rm pri-scc} - \dot{\mathbf{m}}\_{\rm scc} \, \mathbf{c}\_{\rm pri}\right) \mathbf{h}\_{\rm fg} \tag{14}$$

$$\mathbf{E}\_{\mathbf{j}} = \mathbf{h}\_{\mathbf{j}} - \frac{\mathbf{P}}{\mathbf{p}\_{\mathbf{j}}} + \frac{\left(\mathbf{u}\_{\mathbf{j}}^2 + \mathbf{v}\_{\mathbf{j}}^2 + \mathbf{w}\_{\mathbf{j}}^2\right)}{2} \tag{15}$$

In Equation (15), hj can be correlated as hj <sup>=</sup> hsse,j <sup>+</sup> <sup>T</sup> Tref cp,jdT, where hsse,j is defined as the phase standard state enthalpy and hfg is the vaporization latent heat, which is the difference between water vapor and water-sensible heat. The Evaporation–Condensation model is utilized in this work to describe the mass transfer rate between phases. This model relies on the relation of the phase temperature inside the cell to the saturation temperature. Based on this temperature comparison, a mass transfer occurs with releasing or gaining the latent heat in the case of condensation or evaporation, respectively [34]. The mass transfer rate in the case of condensation, as in the present case, is calculated by Equation (16).

$$\dot{m}\_{\rm V-W} = \chi\_{\rm c\rm s} \propto\_{\rm v} \rho\_{\rm v} \left( \frac{T\_{\rm sat} - T\_{\rm cell}}{T\_{\rm sat}} \right) \tag{16}$$

χcs is defined as the condensation parameter and can be calculated by Equation (17). In addition, it is assumed that the condensation from water vapor to water liquid generates a droplet of spherical shape [29,35].

$$\chi\_{\rm cs} = \frac{6}{\mathrm{d\_{drop}}} \beta\_{\rm cs} \sqrt{\frac{\mathrm{M}}{2\pi r \mathrm{T\_{sat}}}} \left[ \frac{\rho\_{\rm w} \mathrm{h\_{fg}}}{\rho\_{\rm w} - \rho\_{\rm v}} \right] \tag{17}$$

In the above equation, ddrop is the diameter of the water droplet and it can be determined by dividing the surface area of the droplet by its volume, as shown in Equation (18) [29,35].

$$\mathbf{A}\_{\text{W}\text{-drop}} = \frac{\pi \mathbf{d}\_{\text{drop}}^2}{\frac{\pi \mathbf{d}\_{\text{drop}}^3}{6}} = \frac{6}{\mathbf{d}\_{\text{drop}}} \tag{18}$$

The concentration of interfacial area is calculated by Equation (19) by assuming equal diameters for all generated water droplets due to condensation [35].

$$A\_{\rm intf-drop} = \frac{\epsilon \alpha\_{\rm V}}{\mathbf{d}\_{\rm drop}} \tag{19}$$

βcs in Equation (17) is the coefficient of accommodation and can be computed by Equation (20) after setting the tunable coefficient of condensation ξcs [29,35].

$$
\beta\_{\rm cs} = \frac{2\xi\_{\rm cs}}{2 - \xi\_{\rm cs}} \tag{20}
$$

### *3.4. Solution Procedure*

In the present simulation, the governing equations are solved with the Evaporation–Condensation model in the three-dimensional domain using the commercial package Ansys Fluent [34]. Additionally, a quick discretization and first-order upwind scheme are used for all conservation equations to save computational time. For more accurate results, the third-order MUSCL discretization scheme is used. The coupled scheme is used for pressure-velocity coupling and the pressure is calculated with the presto scheme, as it is recommended for two-phase complex flow. The coupled algorithm solves the momentum and pressure-based continuity equations together which makes the solution robust and more efficient over the segregated solution scheme. The gravitational acceleration is set equal to <sup>−</sup>9.81 m/s2 in the vertical direction to account for the buoyancy effect with the density gradient. To accomplish robust convergence, various pseudo time values are tested for the fluid and solid zones of the computational domain. The selected pseudo time scale factors for fluid and solid zones are 0.7 and 1, respectively. The pseudo-transient fashion is recommended for the complex two-phase flow, which has high mass transfer through simulation [29]. A grid independence study is performed to grantee the computational results and the generated mesh is shown in Figure 3. The total number of the used tetrahedral mesh is 534,136. The mesh is refined near the tube wall and y+ is approximately equal to 1. Moreover, enhanced wall treatment is used for a smooth transition between the boundary layers. The enhanced wall treatment is recommended in the case of heat transfer simulation and it needs fine mesh in the viscous sublayer.

**Figure 3.** The computational mesh of the inverted U-tube.

Boundary conditions specify the flow and thermal variables on the physical domain boundaries. Figure 2 shows the boundary conditions utilized in the present study. The mass flow inlet applied to each phase and mixture temperature are specified at the flow inlets. The pressure outlet boundary condition is specified as a constant value and is equal to zero-gauge pressure at the pipe outlet. The isolated wall of the duct has been set as hydraulically smooth walls with non-slip boundary conditions.

### **4. Results and Discussion**

In the proposed framework, low energy consumption with a simple construction device for maintaining the low generated pressure for engineering applications is presented. This flow is captured using Ansys Fluent with some model adjustments. The full details for the computational model are shown in Figures 2 and 3. Furthermore, the free predicted water mass flow rate in the U-tube shown in Figure 2 is 3.72 kg/s at a water elevation difference of 75 cm between the two tanks while the mixture inlet at the side pipe is blocked. This free predicted water flow rate is a function of the elevation difference between two water levels. The main objective of the current study is to investigate the two-phase flow interactions when it encounters condensation in the considered U-tube. Besides, the effect of water flow rate ratios and the critical operational conditions are highlighted.

First, the computational results are validated with the results obtained in the experimental part of the present work and with the experimental data, reported in [18] as well, for air and water two-phase flow in pipes, as shown in Figure 4. It is to be noticed here that, in carrying out the computation for validation, the values of the mass flow rate of water siphon flow and air streaming from the side pipe were taken from the present experimental record. Additionally, the air-water two-phase flow in the inverted U-tube is previously studied in [31] and the main findings encouraged the present computations. The air mass fraction in the air/water mixture is represented on the horizontal axis while the air volume fraction appears on the vertical axis. Figure 4 shows an acceptable qualitative validation for the present computational model. It is expected that the water mass flow rate affected by the siphon flow would decrease as the air and water vapor mass flow rate is raised. This is ascribed to the increase in flow resistance. The effect of side air entrainment on the water siphon flow at different water levels between tanks (Δ*Z*) appears in Figure 5, as obtained by the aid of the experimental data in this work. Side air suction from the branch side tube lowers the water siphon flow rate at different water levels between tanks. This is attributed to the effect of entrained air on the vacuum pressure at the tube upper region which is found weakening the siphon flow. The experimentally recorded vacuum pressure at the inverted U-tube tip is shown in Figure 6. Side air entrained into the inverted U-tube decreases the generated vacuum pressure at its tip.

The computational model is executed for the operational conditions of water mass flow rate of 2, 4, 6 and 8 kg/s with side gas mixture fraction by volume of 0.2–0.7. The water vapor fraction in the side mixture is allowed to be in the range of 0.1–0.5. The air mass flow rate effect on the static pressure in the U-tube upper region is shown in Figure 7 at various water vapor mass ratio, MR. Reducing the U-tube upper part static pressure raises the suction effect for the side mixture, which results in improving the present tube function. Increasing both air and water vapor mass flow rates elevates the static pressure in the upper part of the tube, which is undesirable. Additionally, the lowest generated static pressure can be found at a water flow rate of 2 kg/s at lower than or equal to the water/air mass ratio of 1000. For a water/air mass ratio of more than 1000, raising the water vapor flow rates magnifies the static pressure; this is attributed to increasing the water vapor flow rate at a higher air mass flow rate, which decelerates the main flow and pushes the side flow down in the tube and, hence, the condensation takes place away from the upper region. The gradient of increasing the pressure with air mass flow rate in the upper region is more pronounced at a water mass flow rate of 2 kg/s, as noticed in Figure 7a, compared to the static pressure in Figure 7b at a water mass flow rate of 4 kg/s. This is attributed to raising the water mass flow rate, which magnifies the condensation occurrence around the upper region and hence the static pressure decreases due to phase change. The critical flow

rate ratio has the largest value, at which the flow of water caused by the siphon effect withdrawing air and water vapor from the side tube is maintained.

**Figure 4.** Comparison of the CFD results with those of the present experimental ones and reported in [18].

**Figure 5.** Effect of side air on the water siphon flow.

**Figure 6.** Effect of side air on the vacuum pressure at inverted U-tube tip.

**Figure 7.** Effect of air and vapor mass flow rates on the static pressure generation in the tube upper region.

The turbulent intensity and kinetic energy at a point measure the degree of interactions and interruption due to fluids mixing. Figures 8 and 9 record the turbulent values due to mixing the water liquid with the side mixture at the tube upper region. The turbulent intensity and kinetic energy increase with air mass flow rate until . ma = 0.002 kg/s with a water to air mass ratio of 1000 and water mass flow rate of 2 kg/s, as depicted in Figures 8a and 9a. Increasing the water mass flow rate to 4 kg/s has no significant effect on these turbulent values at the tube upper region as, by increasing the water liquid mass flow rate, the side entrained mixture is directed to flow down the tube. So, the generated turbulent values due to the mixing process in Figures 8 and 9 are less noticeable at the tube upper region at water mass flow rates of 4 kg/s compared to 2 kg/s. That means that working at a water siphon mass flow rate of 2 kg/s requires a water/air mass ratio of 1000 or less for avoiding system flow discontinuity. This can be ensured by tracking the flow velocity inside the upper region at the same operational conditions as demonstrated in Figure 10. The flow velocity rises with increasing the air mass flow rate until a mass ratio of 1000 then it recedes, and this inverted tube breaks down water siphon flow at a vapor mass ratio of 0.5, as shown in Figure 10a.

**Figure 8.** Effect of air and vapor mass flow rates on the turbulent intensity in the tube upper region.

**Figure 9.** Effect of air and vapor mass flow rates on the turbulent kinetic energy in the tube upper region.

**Figure 10.** Effect of air and vapor mass flow rates on the flow velocity in the tube upper region.

The effect of flow mass ratios of water, water vapor, and air on the static pressure contours are visualized in Figures 11–14. Low values of air and vapor flow rates satisfy the minimum pressure inside the U-tube, see Figure 11. On the other side, raising the side mixture volume ratio to 0.7 with preserving the vapor mass ratio to be 0.1 increases the static pressure at the top of the U-tube as demonstrated in Figures 11 and 12. Raising the vapor mass ratio from 0.1 to 0.5 with maintaining the side mixture at 0.2 reduces the static pressure generation along the inverted U-tube as depicted in Figures 11 and 13. Meanwhile, increasing both ratios of water and air magnify the static pressure inside the considered tube, but water vapor has the noticed effect on preserving the lowest static pressure as preferred for this tube operation compared to full side mixture effect, as shown in Figures 12 and 14. This means that more condensation, due to increasing the vapor flow rate, enhances the present desired operations of this U-tube and it magnifies the required induced vacuum.

**Figure 11.** Pressure contours at different flow ratios and Xw = 4306.59, α<sup>m</sup> = 0.2 and Xv = 0.1.

**Figure 12.** Pressure contours at different flow ratios and Xw = 461.42, α<sup>m</sup> = 0.7 and Xv = 0.1.

**Figure 13.** Pressure contours at different flow ratios and Xw = 6187.64, α<sup>m</sup> = 0.2 and Xv = 0.5.

**Figure 14.** Pressure contours at different flow ratios and Xw = 662.96, α<sup>m</sup> = 0.7 and Xv = 0.5.

The effect of flow mass ratios in the tube on the side air distribution is shown in Figures 15–18. Increasing the main water flow rate moves the side entrained air to flow downward, not to the suction region in the upper part of the U-tube. This is attributed to force balancing among the inertia of water flow and gravity with the force due to flow suction in the upper portion of the tube. So, reducing the force due to flow suction by magnifying the resulted force due to flow inertia weakens the entrainment operation in the inverted U-tube, as expected above in Figure 9. The flow inertia is a function of the flow velocity, which is determined by the water flow rate. This finding ensures that the relation between the flow rates of water in the tube and water vapor in the side mixture controls the entrainment process of the tube. Moreover, raising the side air flow rate converts the flow in the right branch of the U-tube to bubbly flow which is found to break down the entrainment process and siphon the flow of the U-tube by raising the static pressure in the tube upper part. The water vapor contours are seen in Figures 19–22. Lower values of water flow rates even at higher values of the side water

vapor flow decrease the pressures in the U-tube upper part. This is attributed to the instantaneous local condensation of water vapor in the upper part which enhances the continuous entrainment of the side mixture. Additionally, the water vapor can be seen at the right branch of the U-tube exit without condensation at higher values of water and side mixture flow rates. This process with time breaks the siphon flow in the proposed system.

**Figure 15.** Air volume fraction contours at different flow ratios and Xw = 4306.59, α<sup>m</sup> = 0.2 and Xv = 0.1.

**Figure 16.** Air volume fraction contours at different flow ratios and Xw = 461.42, α<sup>m</sup> = 0.7 and Xv = 0.1.

**Figure 17.** Air volume fraction contours at different flow ratios and Xw = 6187.64, α<sup>m</sup> = 0.2 and Xv = 0.5.

**Figure 18.** Air volume fraction contours at different flow ratios and Xw = 662.96, α<sup>m</sup> = 0.7 and Xv = 0.5.

**Figure 19.** Vapor volume fraction contours at different flow ratios and Xw = 4306.59, α<sup>m</sup> = 0.2 and Xv = 0.1.

**Figure 20.** Vapor volume fraction contours at different flow ratios and Xw = 462.42, α<sup>m</sup> = 0.7 and Xv = 0.1.

**Figure 21.** Vapor volume fraction contours at different flow ratios and Xw = 6187.64, α<sup>m</sup> = 0.2 and Xv = 0.5.

**Figure 22.** Vapor volume fraction contours at different flow ratios and Xw = 662.96, α<sup>m</sup> = 0.7 and Xv = 0.5.

Generally, in mixing two different flows, the turbulent Reynolds number is computed to capture the degree of interaction and turbulence due to mixing. This number is a relation between the turbulent kinetic energies to the dissipation rates and is defined as (*Re*<sup>t</sup> = k2/(εϑ)). Moreover, the turbulent Reynolds number clarifies the mixing shear layers of the flow. The turbulent Reynolds number contours at different values of flow mass ratios are depicted in Figures 23–26. It is noticed that at a water flow rate of 2 and 4 kg/s and side mixture volume ratio of 0.7, higher values of turbulent Reynolds number and flow eddies occur. This ensures, better entrainment, and mixing can be satisfied with the flow conditions of 2–4 kg/s with a mass ratio in the range of 461–662 or less. Acceptable entrainment can work until a mass ratio of 1000 at . mw <sup>=</sup> 2 kg/s, as predicted previously in Figure 10.

**Figure 23.** Turbulent Reynolds number contours at different flow ratios and Xw = 4306.59, α<sup>m</sup> = 0.2 and Xv = 0.1.

**Figure 24.** Turbulent Reynolds number contours at different flow ratios and Xw = 461.42, α<sup>m</sup> = 0.7 and Xv = 0.1.

**Figure 25.** Turbulent Reynolds number contours at different flow ratios and Xw = 6187.64, α<sup>m</sup> = 0.2 and Xv = 0.5.

**Figure 26.** Turbulent Reynolds number contours at different flow ratios and Xw = 662.96, α<sup>m</sup> = 0.7 and Xv = 0.5.

The radial distribution of air and vapor volume fraction contours at different mass ratios and positions (L/d = 10, 25, 50, and 100) are shown in Figures 27 and 28. In these Figures, L represents the measured downward distance on the tube from the side pipe centerline, while d is the tube diameter. Air is flowing near the side mixture inlet at closer positions of the side pipe due to the water steam inertia. On the other side, the vapor is instantaneously condensed after mixing with the tube water stream at lower values of water mass flow rate. Increasing the tube water flow rate makes the vapor continuing downward in the tube without condensation and in some extreme cases, the vapor may exit the inverted U-tube in the lower water tank without condensation. In addition, the vapor is attracted to flow in the tube center cavity compared to air, which is found to flow closer to the tube walls.

**Figure 27.** Radial air volume fraction contours at different downward distances from the side pipe centerline at various flow ratios and Xw = 662.96, α<sup>m</sup> = 0.7 and Xv = 0.5.

**Figure 28.** Radial vapor volume fraction contours at different downward distances from the side pipe centerline at various flow ratios and Xw = 662.96, α<sup>m</sup> = 0.7 and Xv = 0.5.

### **5. Conclusions**

Using the inverted U-tube for side gas entrainment is experimentally and numerically investigated in this paper. Various gases can be entrained from the required cavities, but the air is suggested in this paper as it is the dominant gas in most applications. Side entrainment is considered a mixture of air and water vapor with different mass ratios. The CFD model is validated first with experimental data and an acceptable comparison is obtained. The present experiment shows the effect of mixing side-air with water flowing in the inverted U-tube on the created water siphon flow at various water levels between water inlet tanks due to altering the generated vacuum at the tube tip. Additionally, the siphon water mass flow rate is found to be reduced with increasing side-air mass flow rate. Then, different mass ratios for the main water stream and side mixtures in the inverted U-tube are examined. The proposed U-tube is capable to bring down the absolute static pressure to.55.1 kPa in the tube upper part. This low-pressure region motives the suction of the side mixture to the tube cavity. This study proves that, at water streams of 2–4 kg/s, better side entrainment can be found at the water/side mixture (air and water vapor) mass ratio of 461–662. In addition, the side entrainment can still exist at a mass ratio of 1000 for water flowing at 2 kg/s. Moreover, the condensation can magnify the generated vacuum at the tube upper portion at a higher water mass flow rate. This means that circulating flow in the range of 2–4 kg/s in the U-tube enhances the function of the tube in entraining the side mixture with maintaining the siphon flow and undergoing vapor condensation. Increasing the water mass flow rate weakens the side entrainment of gases to the tube upper part and hence reduces the tube function at certain mass ratios. These findings can be applied to a variety of engineering applications for evacuating systems.

**Author Contributions:** Conceptualization A.H., A.E. and K.Y. formal analysis, K.Y., A.E. and A.H.; funding acquisition A.E.; investigation, A.H., K.Y. and A.E.; methodology, K.Y., A.H. and A.E.; project administration, A.E.; Supervision, A.H. and A.E.; writing—original draft, K.Y.; writing—review and editing, A.H., K.Y. and A.E.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

### **References**


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

© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **In-Vitro Experimental Modeling of Oscillatory Respiratory Flow in a CT-Scanned OSAHS Tract**

### **Zhenshan Zhu 1, Yaping Ju <sup>1</sup> and Chuhua Zhang 1,2,\***


Received: 19 October 2020; Accepted: 9 November 2020; Published: 10 November 2020

**Abstract:** Obstructive sleep apnea-hypopnea syndrome (OSAHS) is a highly prevalent respiratory disorder. The knowledge of respiratory flow is an essential prerequisite for the establishment and development of OSAHS physiology, pathology, and clinical medicine. We made the first in-vitro experimental attempt to measure the oscillatory flow velocity in a computed tomography (CT) scanned extra-thoracic airway (ETA) model with OSAHS by using the particle image velocimetry (PIV) technique. In order to mimic respiration flow, three techniques were adopted to address difficulties in in-vitro experimental modeling: (1) fabricating the obstructive ETA measurement section with the CT-scanned data of an OSAHS patient airway; (2) maintaining the measurement accuracy by using the optical index-matching technique; (3) reproducing the oscillatory respiratory flow rates with the compiled clinical data of transient tidal volumes. The in-vitro measurements of oscillatory respiratory flow velocity manifested the time evolution of the complex OSAHS flow patterns, and the potential wall collapse of the ETA model with OSAHS.

**Keywords:** in-vitro experimental modeling; CT-scanned OSAHS model; oscillatory respiratory flow rates; particle image velocimetry technique; optical index matching

### **1. Introduction**

The human airway is a key organ for the body to exchange air with the external environment. However, this passageway is susceptible to inhaled particulate matters, extreme atmospheric conditions, bad individual habits, and natural aging, resulting in a wide spectrum of respiratory diseases ranging from the simple common cold, snoring, nasal airway blockage and obstructive sleep apnea-hypopnea syndrome (OSAHS) [1]. In the case of such disorders, OSAHS is highly prevalent and poses a high risk of damaging human health and life [2–5]. Its severity is highly correlated to the flow characteristics of the upper airway [6]. Thus, the knowledge of OSAHS respiratory flow in the upper airway is of scientific significance for the establishment and development of OSAHS physiology and pathology, and of clinical medicine references for the diagnosis and therapy of OSAHS diseases.

Compared with the in-vivo measurement of respiratory flow in a realistic human airway, the in-vitro measurement of respiratory flow can effectively provide more detailed flow fields, allowing the respiratory physiologist to deepen the understanding of OSAHS pathology and the respiratory physician to improve diagnosis and therapy for OSAHS diseases. Nonetheless, OSAHS physiological studies of airway anatomy has earned much less attention than healthy airways. Only a small number of studies have focused on respiratory flow in OSAHS-like tract, although much experimental research on the respiratory flow has been carried out in simplified smooth multiple-generation airway models, which were briefly reviewed in Zhu et al. [7]. Following a series of numerical works on

steady inspiratory or expiratory flow in smooth symmetrical and asymmetrical multiple-generation bifurcating airway models [8–10], Liu and his collaborators extended their numerical methods to unsteady respiratory flow in a human upper airway (HUA) model before and after OSAHS surgery [11]. Their numerical results confirmed that after surgery, the HUA flow resistance decreases with even more streamlined flow patterns, demonstrating the potential application of OSAHS respiratory flow in OSAHS surgery. Later, they numerically simulated the steady and unsteady flows in a HUA model with OSAHS resulting from nasopharyngeal obstruction [12]. Their numerical results indicated that the respiratory flow with OSAHS features strong flow injection caused by narrowing of the pharynx at both inspiration and expiration stages. Recently, Song et al. [13] numerically simulated steady inspiration flows in a nasopharyngeal obstruction model with OSAHS. Their model is quite similar to that used in Liu et al. [12]. The numerical results found that inspiratory flow velocity reached a maximum at the obstructed area with a strong negative pressure gradient. Very recently, Cui et al. [14] used the large eddy simulation (LES) method to investigate oscillatory respiration flows in an idealized extra-thoracic (ETA) airway. The simulated results indicated that the unsteady airflow structures are highly impacted by inspiration and expiration phases, and the techniques to model the breathing process.

As for the experimental measurements of respiratory flow, Kim and Chung [15] used tomographic particle image velocimetry (PIV) techniques to measure the airflow in disordered and corrected nasal cavity models. Later, Kim and his collaborators evaluated the feasibility of 3D steady tomographic PIV measurements in a normal nasal cavity [16]. More recently, Wu et al. [17] implemented a smoke-wire visualization of steady flow at the mid-sagittal plane of a HUA model with OSAHS by using a high-speed camera. Their experimental visualization, together with numerical simulation, suggested that for OSAHS patients, HUA narrowing is a symptom of OSAHS development and deterioration.

All of the above in-vitro experimental models for OSAHS respiratory flow were limited to steady flow cases. The recent numerical observations [18,19] of the unsteady flow patterns differ significantly from the steady ones in the ETA models. In-vitro measurement of unsteady respiratory flow in human airway models poses a great challenge for fluid experimentalists due to related geometrical complexity, limited optical accessibility and oscillatory flow rate [20]. Specifically, three difficulties need to be overcome if an accurate and reliable mimicking of respiratory flow in in-vitro environment was required, i.e., a realistic respiratory tract model with OSAHS, good optical index matching, and reliable control over the oscillatory respiratory flow rate. It is worth noting that most investigations on oscillatory respiratory flow were based on the simplified sin curve flow rate as the oscillatory boundary condition [18,19]. In fact, the realistic oscillatory respiratory boundary condition is quite different from the sin function [21].

The goal of the present research is twofold. The primary goal is to tackle the above mentioned three critical difficulties arising in the in-vitro experimental modeling of OSAHS flows. The secondary goal is to employ the experimental rig to measure the oscillatory respiratory flows in an ETA model with OSAHS. The ETA model with OSAHS is composed of an oral cavity, pharynx, larynx and trachea with total closure of the nasal passageway.

The remainder of this paper is organized as follows. In Section 2, an in-vitro experimental model will be described, focusing on techniques to overcome the above three difficulties. Section 3 will present PIV techniques tailored to the experimental rig. Section 4 will present the PIV measurement results of oscillatory respiratory flow in the ETA model with OSAHS and their clinical implications. Concluding remarks will be made in Section 5.

### **2. In-Vitro Experimental OSAHS Flow Model**

### *2.1. ETA Model with OSAHS*

Creating an anatomically correct transparent experimental respiratory flow model is essential to measure respiration flow inside a human airway by using the PIV technique. The present ETA model was based an anatomical data, which was obtained from a scanned airway of a 48-year-old Chinese male OSAHS patient with a 128-slice spiral CT, under the rest breathing condition. The apnea-hypopnea index (AHI) and body mass index (BMI) were 52 and 26.8, respectively. In order to obtain detailed data, CT scanning was performed on the patient along the coronal, sagittal, and cross-sectional directions. The thickness between the CT-scanned slice images was 1 mm. The OSAHS condition of the patient was so severe that the passageway through the nasal cavity was totally blocked in the retro uvula region (RUL) as shown in Figure 1. The ETA model composed of the oral cavity, pharynx, larynx and trachea was adopted to mimic the patient's breathing through the oral cavity. The ETA model with OSAHS was then constructed with the image processing software Mimics [15,16,22]. The digital ETA model is illustrated in Figure 2, from which the obstruction resulting from the retro uvula hypertrophy can be clearly observed in Figure 2a. Two representative measurement sections were selected, i.e., the global A-A section viewed from the lateral side and the local B-B section near the obstructed region, as shown in Figure 2a,b.

**Figure 1.** Sagittal plane CT image of OSAHS patient.

**Figure 2.** Reconstructed ETA model with OSAHS: (**a**) dorsal view; (**b**) lateral view; (**c**) front view; (**d**) four cross-sections.

To facilitate an accurate experimental measurement, the experimental model was finely fabricated as follows. First, a negative ETA model was printed with water-soluble material using a fused deposition modeling (FDM) printer [16,23,24]. The melting temperature of the water-soluble material was 220 ◦C and the wire diameter was 1.75 mm. The maximum printing error of the negative model was less than 0.1 mm. Second, the negative model was coated with a water-soluble glue. The coating glue functioned as a slice separator of the negative model as the glue dried out. Both the water-soluble material and glue were mainly composed of polyvinyl alcohol (PVA), but they had different physical properties. The dried water-soluble glue will not stick to the external filler silicone, but the water-soluble material will. Third, the coated negative model was put into an optical glass box that was then carefully filled with transparent liquid silicone Sylgard 184 [25]. Before it was poured into the box, the silicone was needed to be vacuumed to remove bubbles in it. If not, the fabricated model wall would absorb tiny bubbles, which then seriously polluted measurement results. The liquid silicone was solidified by curing at room temperature over 48 h. Finally, the negative model and the coating glue were washed with cold water. The fabricated experimental model is shown in Figure 3. The compliance of the human airway tract and the mucus attached to the inner wall [26,27] were not considered in the present experimental model for the sake of accurate PIV measurements.

**Figure 3.** Optical index matching between ETA model and working fluids: (**a**) air; (**b**) water; (**c**) mixture of water and glycerol.

### *2.2. Optical Index Matching*

In optical measurements of the flow inside the airway model, the optical index matching between the working fluid and model wall is crucial [20,28–31]. Good index matching can eliminate optical distortion to ensure measurement accuracy. After a series of priori tests with different mixtures of water and glycerol as working fluid, the mixture composed of 60.7% glycerol and 39.3% water [32] reached the best optical index matching. For comparison, Figure 3a–c presents optical index matching with air, water and water-glycerol mixture, respectively, as working fluids in the model. Clearly, the best optical index match was achieved with the water-glycerol mixture. As shown in Figure 3c, the lines in the quadrille paper sheet do not suffer from deformation as one looks through the model [20].

### *2.3. Oscillatory Respiratory Flow Rate*

In the present unsteady measurement, the realistic normal male breathing flow rate transients at rest were used to determine the modeling flow rate transients. Statistical data of the breathing flow rate transients compiled by a Beijing medical institution from hundreds of supine male adults were used. Similar flow rate transient curves were also confirmed in Kim and Chung [21] and Huang and Zhang [33]. The tide volume (*V*T) here was chosen to be 500 mL. The case where the inspiratory air volume is a little larger than the expiratory one [34] was not considered in the present test rig. The typical breathing cycle period at rest is within 3–5 s; 5 s were chosen due to the severe status of the OSAHS patient.

To mimic the realistic flow physics and respiratory physiology, the similarities of Womersley number (α = *L*/2· 2π· *f*b/ν) and Reynolds number (*Re* = *V*·*L*/ν) were kept between the breathing air and the modeling water-glycerol liquid [23,29,32,35]. The two similarity criteria ensured flow similarity such as the same relationship of pressure drop coefficient vs. flow coefficient. The Womersley number is a dimensionless number in biofluid mechanics representing the ratio of transient inertial force to viscous force and is an important similarity criterion for the oscillatory respiratory flows. Here, the *Re* number can be defined as *Re* = 4*f* <sup>b</sup>·*V*T/(*L*·ν), as mentioned by Janke et al. [23,24]. In these definitions, *L*, *V*T, *f* <sup>b</sup> and ν refer to the inlet hydraulic diameter of the tract, tide volume, frequency of the oscillatory respiration, and kinematic viscosity of the working fluid, respectively. The kinematic viscosity (ν) of the water-glycerol mixture and breathing air was 8.933 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>m</sup>2/s and 15.03 <sup>×</sup> <sup>10</sup>−<sup>6</sup> m2/s, respectively, at room temperature during the experiments. In this oscillatory respiratory flow, the maximum Reynolds numbers during inspiration and expiration phases were 1442 and 1149, respectively, and the Womersley number was 2.89. To maintain the similarity of Reynolds number and Womersley number, the breathing cycle period (1/*f* b) of 5 s was tuned to be 8.5 s in the present in-vitro experimental model where the working fluid was taken as the water-glycerol liquid. The breathing flow rate transients as well as those employed in the experiment are shown in Figure 4. Under such oscillatory boundary condition, roles of inlet and outlet of this model would switch periodically. At the inspiration stage, the boundary at the mouth side can be considered as the inlet, and the boundary at the trachea side was the outlet. While at the expiration stage, the contrary would be the truth.

**Figure 4.** Human breathing and modelling flow rate transients.

To mimic the modeling flow rate transients as shown in Figure 4, several apparatuses and procedures were carefully designed in this research. The schematic and photo of the experimental setup are presented in Figure 5. A programmable logic controller (PLC), servo-motor, servo-motor driver, gear box, bent axle, piston, air exhausting device, reservoir and *V*<sup>T</sup> monitor (Figure 5a) were employed to provide precise control over piston movements and thus experimental flow rate transients. Figure 5b,c presents a photo of the main apparatus and experimental rig, respectively. First, a PLC program was developed and run on a personal computer to transfer the PLC input digital signal into the PLC output analog voltage signal. Second, step by step, the driver of the servo-motor drove the movements of the servo-motor. The duration of the movement in each step was 20 ms. Thus, a total of 425 step movements were precisely controlled over the present single respiratory cycle (8.5 s). Third, the gear box was used to increase the servo-motor torque output, which eventually actuated the bent axle and the piston. Lastly, the air exhausting device was used to remove potential air bubbles absorbed into the test rig. The apparatus of the reservoir and *V*<sup>T</sup> monitor was mainly made

of a graduated cylinder, as shown in the Figure 5b. This cylinder can not only collect working fluids and add tracer particles, but also monitor the output *V*T. The liquid level in the graduated cylinder would fluctuate up and down, following simultaneously with the oscillatory experimental boundary. By recording the locations of the highest and lowest liquid levels in the graduated cylinder, the output *V*<sup>T</sup> was monitored.

(**a**)

(**b**)

**Figure 5.** Experimental rig: (**a**) schematic; (**b**) main apparatus (except ETA tract and PIV system):1 PLC; 2 driver; 3 servo-motor; 4 gear box; 5 bent axle; 6 piston; 7 air exhausting device; 8 reservoir and *V*<sup>T</sup> monitor; (**c**) photo of experimental rig.

(**c**)

### **3. PIV Techniques for the ETA Model with OSAHS**

The other goal of this research was to apply the above experimental rig to measure the oscillatory respiratory flow velocity in an ETA model with OSAHS. The LaVision 2D-2C (two-dimensional, two-component) PIV system was used to measure the unsteady flow velocity at the two representative sections, i.e., section A-A and section B-B, as shown in Figure 2a,b.

A double pulse yttrium-aluminum-garnet (YAG) laser and a high-speed charge-coupled device (CCD) camera were employed to obtain digital images. A laser (Evergreen Corp.) with power up to 600 mJ per pulse was used to provide double laser sheets with frequency of 15 Hz, wavelength of 532 nm and thickness of approximately 1 mm to induce fluorescence of tracers. The CCD camera (Imager SX 4M) was used to store digital images. The camera features an interline transfer chip with progressive-scan low-noise readout, high spatial resolution of 2360 × 1776 pixels (4 Mpixel, 5.5 <sup>×</sup> 5.5 <sup>μ</sup>m<sup>2</sup> per pixel size) in 12-bit mode, and fast frame rates up to 31 Hz. Here, 15 Hz frame rates were used to match the frequency of the laser pulse.

Determining the distinctive time (*dt*) delay of two subsequent particle images (double shutter) in the PIV measurements of oscillatory flow is an important task since the *dt* setting must be adapted to the maximum flow speed as well as the flow reversal condition [32]. Here, *dt* was determined according to the CCD spatial resolution, field of view of the camera, estimated average velocity, and target pixel shift of tracing particles in double shutter, which were 2360 <sup>×</sup> 1776 pixels, 15 <sup>×</sup> 11.3 cm2, 0.7 m/s, and 7 pixels, respectively. The *dt* and the duration of exposure time of particle images were determined to be 500 μs and 200 μs, respectively. The length of camera image (11.3 cm) was smaller than that in ETA model (approximately 20 cm); consequently, one PIV measurement of A-A section had to be divided into two sub-sections. The two digital images obtained from the two sub-sections measurements at the same phase were then stitched to form a whole PIV measurement image. Within the present modeling cycle (8.5 s), 128 groups of double shutter images were continuously acquired with the PIV system and controlled with a synchronizer. All steps of control, acquisition and post processing were carried out in a DaVis 10.0.5 software environment [36].

The reflected light is very harmful to both the camera CCD itself and the CCD acquisition quality of near-wall flows. To avoid light reflection from walls, good optical index-matching, suitable tracer fluorescence and a long-pass glass filter were employed [16]. The working fluid with good optical index-matching not only eliminates optical distortions, but also reduces wall light reflection. Fluorescent tracers were used to scatter 590 nm wavelength light under 532 nm wavelength laser irradiation. Using the long-pass glass filter before the camera, the reflected light with wavelength below 550 nm will be filtered out, thus greatly alleviating wall reflection.

The long-pass glass filter and fluorescent tracers used in PIV measurement will not distort measured results. The long-pass glass filter had more than 99% transmittance of the scattered light with wavelength above 550 nm. The mean diameter *d*<sup>p</sup> and density ρ<sup>p</sup> of fluorescent tracers were 4 μm and 3.4 g/cm3, respectively. The response time (τp) of tracing particles defined as *d*<sup>p</sup> <sup>2</sup>ρp/(18μ) was calculated as 2.93 <sup>×</sup> <sup>10</sup>−<sup>7</sup> s [30]. The dynamic viscosity <sup>μ</sup> of the water-glycerol mixture was 1.03 <sup>×</sup> <sup>10</sup>−<sup>2</sup> Pa·s.

### **4. Results and Discussion**

In a respiration cycle, four representative phases were highlighted to present PIV measurements of oscillatory flow velocity fields at two typical A-A and B-B sections, as shown in Figure 6. Locations of A-A and B-B sections are defined in Figure 2a,b. The four selected phases represent the reversal phase of expiration-inspiration (phase one), the peak phase of inspiration (phase two), the reversal phase of inspiration-expiration (phase three), and the peak phase of expiration (phase four), respectively. The measured flow velocity contours were colored with the magnitude of absolute velocity (ABS velocity). The measured results will be presented and discussed in the following two sub-sections for reversal phases (phases one and three) and peak phases (phases two and four). For a clear presentation, local views at B-B sections were deliberately enlarged. In order to quantitatively analyze measured results, the main-stream velocity transients at four cross-sections were also presented. The locations and the main-stream velocities profiles of the four cross-sections are shown in Figures 2d and 7, respectively.

*Appl. Sci.* **2020**, *10*, 7979

**Figure 6.** *Cont.*

**Figure 6.** Measured results at four representative phases at A-A section and B-B section: (**a**) phase one: expiration-inspiration reversal transient; (**b**) phase three: inspiration-expiration reversal transient; (**c**) phase two: peak inspiration transient; (**d**) phase four: peak expiration transient.

**Figure 7.** Measured main-stream velocities at four cross-sections at four representative phases: (**a**) cross-section 1; (**b**) cross-section 2; (**c**) cross-section 3; (**d**) cross-section 4. Locations of these four cross-sections are defined in Figure 2d.

### *4.1. Reversal Phases*

At phase one, the respiratory flows are undergoing transition from expiration to inspiration (Figure 6a). The flow patterns presented are irregular both at the A-A section and the B-B section. The reversal flow and jet flow structures can be clearly observed at the A-A section. The interaction between the expiration and inspiration flows leads to strong flow reversal at the oral cavity and oropharynx region, which is hardly considered as either expiration or inspiration pattern. However, the local flow structure at the B-B section mainly corresponds to the expiration pattern with a large portion of expiration flow at the mainstream region, except only a small portion of reversal flow at the near wall region. The flow pattern at the expiration-inspiration phase is more complex than that observed in Adler and Brücker [22], who claimed that the flow at this phase could be regarded as expiration flow. This difference may be caused by more realistic oscillatory flow rate boundary conditions used in this research rather than the simplified sine-law flow rate boundary conditions used in Adler and Brücker [22]. Such jet-reversal flow structure is an indicator of complex

expiration-inspiration transition at the quasi-static conditions where the net respiratory flow rate is zero while the unsteady inertial force and thus the respiratory flow acceleration and velocity are not zero.

We can also see that compared with the quasi-static condition at phase one, the jet-reversal flow structure occurs even more obviously at phase three (Figure 6b), due to the milder transition from inspiration to expiration (Figure 4). Neither an inspiration nor an expiration pattern is the dominating flow pattern at phase three. The jet-reversal flow structure extends into the oral cavity, oropharynx, laryngopharynx and trachea (Figure 6b, A-A section), and the widely spread reversal flow vortexes at the obstructed region (Figure 6b, B-B section). The inspiration-expiration pattern is very different from that in Adler and Brücker [22], who considered the flow at this phase to can be regarded as inspiration flow.

From the reversal transients, we could conclude that the realistic oscillatory respiratory boundary condition and the realistic model have a significant effect on respiration pattern, and associated flow structures. This observation can be clearly confirmed in Figure 7. Although the main-stream velocity profiles at reversal phases one and three are relatively small, their patterns are not similar to those at inspiration phase two and expiration phase four. Although quite few numerical and experimental studies focused on the reversal transients, the present PIV measurements confirm that the respiration patterns at the reversal phases can neither be regarded as the pattern of inspiration nor that of expiration due to the occurrence of jet-reversal flow and quasi-static condition.

### *4.2. Peak Phases*

At peak inspiratory transient (phase two), the inspiration flows form an obvious jet-wake structure and a strong flow reversal (Figures 6c and 7b). Jet-wake flow structure, which often occurs in the human trachea and bronchi due to wall curvature [7], also appears at this phase in the ETA model with OSAHS. High-speed jet fluids accumulate at the dorsal side of the laryngopharynx region, while the anterior side of the trachea (Figure 6c, A-A section, Figure 7a,b). The narrowed tract does have a significant impact on inspiration flow at this phase with a clear injection phenomenon. As shown in the enlarged vector diagrams at the phase of peak inspiration, the jet-wake structure together with the flow reversal enhances fluid mixing and heat transfer, which benefits the velocity, temperature and humidity conditioning of the inhaled air as confirmed by Kim et al. [37].

At the peak expiration phase (phase four), however, the expiration flow shows very different flow patterns (Figure 6d). The high speed jet fluids accumulate at the anterior side of the laryngopharyngeal region and the dorsal side of the trachea (Figure 6d, A-A section, Figure 7a,b), which are at opposite sides of the tract at the peak inspiration phase (Figure 6c, A-A section). Such differences may be due to the combined effect of convection and viscous motions in the curved tract and complex wall attachments. The injection phenomenon caused by the narrowed region resembles the numerical observation in Liu et al. [12]. It should be noted that we used the oropharyngeal model in our research rather than the nasopharyngeal model used in Liu's study. The flow patterns at the trachea region are very different with the results in Wang et al. [38] using a simplified model and the results in Cui et al. [14] using an idealized model. In those studies, the expiration flow patterns at the trachea region were quite regular and presented straight tube-like flow. The differences may be caused by the more realistic model used in this research.

From the results above, we can clearly conclude that the inspiration and expiration flow patterns are quite different from each other. High speed jet fluids appeared at the opposite sides of the oropharynx, laryngopharynx and trachea regions, which was not highlighted in the previous similar numerical studies [6,14,19,26,37]. Expiration process is not a simple reversal process of inspiration.

Considering that the patient's nasal passageway was completely blocked, we believe that the collapse of the passageway most likely occurs at the peak inspiration phase at the obstructed region, thus resulting in apnea. There are two main reasons for this: firstly, the effect of obstructed tract on oscillatory respiratory flow fields at peak phases is much more apparent than that at the reversal phases. The occurrence of injection phenomenon at the narrowing obstructed area at the peak inspiration

phase indicates very large pressure differences over this short, obstructed distance [13]. Secondly, at inspiration phase, the driving force induced by the movements of the thoracic cavity and abdominal cavity is a negative pressure. The combination of large pressure difference across the obstructed area and negative pressure at the peak inspiration phase would more easily induce collapse of passageway wall and eventual blockage of the tract, thus resulting in apnea. This condition worsens and becomes more dangerous for the deep-sleeping OSAHS patient, whose passageway tissue walls become very loose.

Overall, we presented a temporal and spatial evolution of the oscillatory respiratory flow patterns within a CT-scanned model with OSAHS using PIV measurements. It should be recognized that the calibrated CFD (computational fluid dynamics) method can provide even more flow details that are hardly, if not impossible, captured with PIV techniques [39], such as pressure drop and airway resistance in the CT-scanned complex geometry. Airway compliance is another important issue in human tract [40] that was not modeled in the present study. However, the present work contributes to an experimental modeling of the complex oscillatory respiratory flow by using PIV measurements and can be regarded as a step-stone for an even deeper understanding of OSAHS airway flows.

### **5. Conclusions**

In the present research, particular efforts were made to experimentally mimic the respiratory flow conditions in an in-vitro experimental model to measure the oscillatory respiratory flow velocity inside a CT-scanned ETA model with OSAHS by using PIV technique. Four representative phases were highlighted to clarify the respiratory flow patterns and the associated flow structures in the ETA model.

With realistic non-sinusoidal respiratory flow rate and CT-scanned ETA model, the respiratory flow patterns at the reversal phases can neither be regarded as the pattern of inspiration nor that of expiration.

Flow patterns at the inspiratory and expiratory phases are very different. High-speed jet flows generally appear at the opposite sides of the oropharyngeal, laryngopharyngeal and tracheal regions at the inspiration and expiration phases. The expiration process is not a simple reversal process of inspiration.

The respiratory flow patterns and the associated flow structures have clinical implications for breathing tract blockage. The area of stenosis is apt to be blocked at the peak inspiration phase due to potential wall collapse.

To the best of the authors' knowledge, this has been the first research devoted to measuring oscillatory respiratory flows inside CT-scanned ETA model with OSAHS. The PIV measurements of unsteady flow velocity will be used to calibrate the computational fluid dynamics method, which can be further used to capture even more detailed flow fields.

**Author Contributions:** Funding acquisition, C.Z., Y.J.; project administration, Y.J.; formal analysis, Z.Z.; investigation, Z.Z.; data curation, Z.Z.; writing—original draft preparation, Z.Z.; writing—review and editing, Y.J. and C.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was financially supported by the National Key Research and Development Project of China under Grant No. 2016YFB0200901, National Natural Science Foundation of China under Grant No. 51176146, and Shaanxi Key Research and Development Project under Grant No. 2018KWZ-01.

**Conflicts of Interest:** The authors declare that for this study they have no conflicts of interest.

### **References**


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

© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Experimental and Numerical Study on Performance Enhancement by Modifying the Flow Channel in the Mechanical Chamber Room of a Home Refrigerator**

### **Dong Kyun Kim**

Division of Mechanical Engineering, Tongmyong University, 428 Sinseon-ro, Nam-gu, Busan 48520, Korea; kimdk@tu.ac.kr; Tel.: +82-10-7402-0744

Received: 6 July 2020; Accepted: 8 September 2020; Published: 9 September 2020

### **Featured Application: Design of high performance refrigerator.**

**Abstract:** Recently, many refrigerators of high have been designed so as to improve the performance of the refrigerator and to economize on energy. There are many methods of improving the total efficiency of the refrigerators using a methodological approach. In this experimental and numerical study, three methods—namely, modifying the flow channel and the operating conditions and changing the positions of the fan at the mechanical chamber room (MCR)—are used to analyze the performance improvement of the refrigerator. These methods of changing the flow channel in an MCR of a refrigerator are, in many ways, used as a regulation of the upper flow region of the condenser, especially in this research. Modifying the refrigerators' MCR is carried out by changing the shape of the cover back machine (CBM). We compute the fluid flow in a refrigerator MCR by computational fluid dynamics (CFD) techniques. The commercial code ANSYS CFX 16.2 is used for this computation. The regulation flow region shows a reduction in power consumption of about 0.75% per month. The result of the operating conditions of the fan at MCR is relative to the flow rate for reducing power consumption. The changing positions are also relative to the length of the reducing power consumption. Moreover, with the results obtained by CFD, we understand the flow structure in a refrigerators' MCR for various types of CBM. The results also show that the efficiency of the refrigerator is improved by 1.2%.

**Keywords:** mechanical chamber room; performance enhancement; flow channel; cover back machine; domestic refrigerator

### **1. Introduction**

Some of the common household appliances are refrigerators, washing machines, clothes dryers, electric ovens, etc. [1]. Among the household appliances, the refrigerator is the most common one due to its indispensable nature and hence, is the most manufactured [1,2]. Currently, there are over a billion domestic refrigerators worldwide [2]. The energy consumption of household appliances increases in proportion to its usage. Among these appliances, the refrigerator is one of the most energy-consuming devices [2]. It is estimated that more than 45% of the world food production would go waste if not for the refrigerators' cold storage and distribution [2]. Thus, improving their energy efficiency, efficient energy consumption, and prolonged life becomes imperative [1,2]. There are different mechanisms to increase the efficiency of the refrigerator, namely, by improving the compressor performance, thermal isolation, heat exchanger design, implementation of optimum control in refrigerator operation [2], improvements in cabinets (such as insulation, improved gaskets), improvement in refrigeration systems (such as fan motors, high efficiency compressors), improved heat exchangers, advanced defrost mechanisms, etc. [1].

Environmental problems, such as global warming, have precipitated a large number of studies—on various factors—that affect the power consumption of a refrigerator to improve their efficiency [3]. On the evaluation of the efficiency and performance of refrigerators, Dmitriyev and Pisarenko investigated the change in the performance of a refrigeration system using a Phase-change Material (PCM) cooling agent [4]. Laguerre et al. investigated the change in the temperature of freezing and refrigeration chambers in a food storage refrigeration system according to its storage temperature and usage frequency in an empirical manner [5]. Sergio Ledesma and J M Belman-Flores mathematically studied the thermal effect of the fresh food compartment with a change in glass shelve locations in a domestic refrigerator. The moving variance and the finite difference approximation were used to analyze the thermal distribution and the temperature change rate, respectively, in their study. The results showed that the use of thermal maps to analyze the mean and variance of the temperature was an important tool to assess a good combination of shelf locations [6]. M Rasti et al. experimentally studied the energy efficiency enhancement of a domestic refrigerator by replacing R134a with R436A without any mechanical alteration. The results showed that by replacing the R436A, the On time and energy consumption was reduced by 13% and 5.3%, respectively, for a day, and the optimum charge for R436A was reduced to 48% when compared to the charge of R134a [7]. Shikalgar and Sapali experimentally studied the energy and exergy analysis of a domestic refrigerator by comparing the hot-wall air-cooled (HWAC) and box type shell and tube water-cooled (BCTWC) condensers. The results showed that the coefficient of performance and the exergy efficiency of BCTWC was increased by 18–20% and 6.89–9.13%, respectively, when compared to the HWAC. In addition, the energy consumption of the BCTWC was reduced by 17% in comparison to HWAC [8]. Belman-Flores et al. experimentally and numerically studied the temperature stratification of the fresh food and freezer compartments of the domestic refrigerator with bottom mount configuration. The results showed that the proposed design had a good distribution of airflow and hence, uniform temperature at the fresh food compartment. Furthermore, in comparison to the original design, the new design showed less ON stages for the compressor [9]. Antonio and Afonso conducted a comparison study between the computational fluid dynamics (CFD) and Artificial Neural Network (ANN) modeling of the air temperature fields inside the refrigerator cabins. The results showed that the ANN model produced a low absolute error of 0.8 K in comparison to 1 K of the CFD model. Moreover, the temperature field obtained with the ANN model was more refined in comparison to the CFD model [10]. A review of advances in domestic refrigerators can be had from [1,2,11]. However, there are very few studies on the improvement in refrigeration performance according to the internal structure of a refrigeration chamber and its conditions that affect the performance of a refrigerator [12–15].

The purpose of this study is to improve the total efficiency of the refrigerator by maintaining the temperature inside the mechanical chamber room (MCR) at a constant level or by decreasing it—based on the improvement of the flow line and the operating condition of the cooling fans—experimentally. In addition, to optimize the flow line inside the MCR by changing the shield shape of the cover back machine (CBM) based on computational flow analysis.

### **2. Factors A**ff**ecting the E**ffi**ciency of the Refrigerator**

The refrigerator consists of a front and rear (F and R) MCR. This chamber is called as a cover back machine (CBM). The MCR is installed at the rear of the refrigerator. Figure 1 shows the component arrangement of the MCR in the refrigerator. It consisted of a compressor, condenser, suction pipe, discharge pipe, cooling fan, guide fan, fan motor, capillary tube, a dryer. The suction pipes were installed at the entrance of the compressor, where a cooling agent moved. The discharge pipes were installed at the exit of the compressor.

In general, the F and R chambers were maintained at an inner temperature of −18 ◦C and +3 ◦C, respectively, where and they showed differences in the outer temperature as 48 ◦C and 27 ◦C, respectively, for an ambient temperature of +30 ◦C [16]. As a normal operation in a refrigerator, the compressor is operated as an On/Off control, and the cooling fan is operated at the same time. A cover back machine (CBM) is installed to prevent the inflow of mechanical noise and other foreign objects into the internal space. The cooling fan is used to cool down the condenser. The cooling fan causes an internal flow inside the refrigerator.

**Figure 1.** Component arrangement of the mechanical chamber room (MCR).

The operational losses, additional device (parasitic) losses, load losses in insulation used to isolate the inflow of external air, input losses in the On/Off control cycles, thermal load losses in the fan motors, input losses in the On/Off control of a compressor, heat exchange efficiency in a condenser, pressure decrease in the pipeline, loss in an evaporator for removing frosts are some of the factors that lead to the inefficiency of a refrigerator [16–19]. In addition, within the MCR, the heat source generated from the compressor, condenser, fans, and motors influence the decrease in the efficiency of the whole cycle of the refrigerator. Moreover, if the heat release is less than the heat released from within the MCR, the temperature of the system will be increased, leading to a decrease in the efficiency of the refrigerator. Thus, a cooling process is applied to the surface of the condenser using cooling fans, and the flow line in the MCR is rearranged to decrease the temperature. However, it is necessary to improve such a flow line to increase the efficiency as there are lots of resistant factors in the chamber due to the complicated flow line.

The internal flow inside the refrigerator caused by a cooling fan is discharged to each vent of the CBM through the forced convection generated by the fan as the compressor condenses refrigerant vapor. The performance of heat transfer inside the MCR caused by such forced convection affects its co-efficient of performance (COP). The COP is the refrigeration effect of a compression work, and it can be calculated as follows [4,14].

$$\text{COP} = \frac{\text{Cooking Efficiency}}{\text{Compressor work}} = \frac{\text{Q}\_{\text{L}}}{\text{AW}} \tag{1}$$

As the vent temperature of a compressor and the condensation temperature of a condenser are decreased at a p-h diagram to improve the COP [18–20], the refrigerant condensed by the condenser may improve the refrigeration effect of COP due to the pressure reduction by an expansion valve. Furthermore, in a convectional heat transfer, by cooling the fan, the heat transfer coefficient of h depends on the internal temperature of the MCR and the air temperature around the vent of the CBM near the cooling fan. Because an effective shield in the vent of the CBM at the front of each major component significantly affects the COP, the convectional heat transfer coefficient, and power consumption, it becomes an important factor in the flow inside the MCR and its flow line [3,13].

#### **3. Background Theory**

The Carnot cycle applied to a refrigerator consists of an isentropic process that corresponds to the work of a compressor, static pressure heat release process in a heat exchanger, expansion process in a capillary tube, and static pressure heat suction process in an evaporator [21]. In general, the static pressure heat release and heat suction process are reversible processes that show a decrease in its efficiency as the temperature difference between high and low temperature regions in a heat exchange process. In addition, it plays a role in decreasing the refrigeration capacity due to the increase in the compression rate of the compressor and the decrease in volume efficiency.

Regarding the solution for solving these problems, although an expansion for the section of heat exchanges is to be recommended, it represents significant difficulties in the structure of the MCR and the space for installing this chamber for the present refrigerator [22]. Therefore, it is the most ideal way to improve the efficiency of a refrigerator by decreasing the temperature in an heat exchanger through changing the flow line using cooling fans, increasing the surface heat transfer coefficient, and decreasing the compression rate of a compressor.

Figure 2 shows the schematic of the reverse Carnot cycle. In the figure, h is enthalpy; h1 and h2 are the initial and final enthalpy state of the compressor, respectively. The numbers 1, 2, 3, and 4 denote the initial enthalpy state of the compressor. The numbers 1- , 2- , 3- , and 4 denote the improved enthalpy state of the compressor. h2 is the improved enthalpy state of the compressor. P is the pressure. AW (= h2 − h1) is the difference between the two enthalpy states. QH and QL are the high and low heat flux, respectively. Thirty degrees denotes the low temperature of the system, and 40 ◦C; 80 ◦C denotes the high temperature of the system. In the operating condition, as shown in Figure 2, the amount of heat transfer and COP can be obtained using the following equations:

$$\text{COP} = \frac{\text{Cooking Efficiency of Evaporator}}{\text{Power Consumption of Compressor}} = \frac{\text{Q}\_{\text{L}}}{\text{W}'} \tag{2}$$

$$\mathbf{Q}\_{\overline{\mathbb{C}}} = \mathbf{K} \mathbf{A} \boldsymbol{\Delta} \mathbf{T} = \mathbf{m}\_{\text{ref}} \boldsymbol{\Delta} \mathbf{h},\tag{3}$$

$$\mathbf{Q\_L = h\_2 - h\_1 \implies h\_{2'} - h\_{1'}},\tag{4}$$

**Figure 2.** Schematic diagram for the reverse Carnot cycle.

The equations above are used to investigate changes in the efficiency of the whole system through improving the heat transfer performance by reducing the temperature difference between high and low temperature regions [23].

### **4. Experimental Conditions**

Figure 3 shows the arrangement of the experimental apparatus used in this study. The experiment was carried out in a constant temperature and humidity chamber to maintain uniform conditions. Two specimens, #A and #B, were used under the same room conditions to ensure the accuracy of the

experiment. In addition, the system consisted of a logger that configures 60 channels for measuring the temperature inside the MCR and the surface temperature of the compressor and condenser, a power meter that measures the power of these specimens, and a static voltmeter that was able to provide a uniform voltage [21].

**Figure 3.** Arrangement of experimental apparatus.

The specification of the constant temperature and humidity chamber and instruments are presented in Table 1. In the experimental conditions, the temperature inside the constant temperature, and humidity chamber and relative humidity were configured as a specific range based on the standard for the test of a refrigerator. Furthermore, a T-Type thermocouple was connected to the logger, the surfaces of the compressor and condenser installed inside the refrigeration chamber and MCR, and inside the MCR to measure the temperature [22,23].


**Table 1.** Specifications of the apparatus.

In addition, the power that was applied to each specimen was measured using a voltmeter in which the output port of the voltmeter was directly connected to the logger. Figure 4 shows the compressor On/Off cycle of the defrost state of the refrigerator. From the figure, it can be seen that the defrost peak occurred twice. The first peak occurred initially at 1 h, and the second peak occurred during the 15th h. The figure also shows the temperature of the refrigerator ranging from 0 ◦C to 20 ◦C. As shown in Figure 4, the internal temperature of the MCR and the surface temperature and power of the compressor and condenser were measured by configuring a stable region between the processes for removing frost. The measurement was repeated five times for each condition and determined as an ensemble average.

**Figure 4.** Compressor On/Off cycle for defrost state.

#### **5. Numerical Analysis**

### *5.1. Governing Equations and Numerical Solutions*

This study was performed under the assumption that the fluid flow type inside the MCR represents a three-dimensional steady state, incompressible, and turbulent flow. The physical properties of the working fluid were kept at a constant level for temperature and pressure for numerical calculation. Due to the disturbance of the fluid by the cooling fan, the flow in the machine room will behave in a complex manner. Therefore, the standard *k* − ε turbulence model that represents excellent estimation performance in an internal flow analysis that has been widely adopted was used to perform flow analysis in the refrigerator machine room [9,10,24]. The governing equations were composed of the continuous equation, the momentum equation, the turbulent kinetic energy equation, and the turbulent kinetic energy dissipation rate equation, and are as shown in the following Equations (5)–(10) [25]:

$$\frac{\partial u\_i}{\partial \mathbf{x}\_i} = 0,\tag{5}$$

$$
\rho \frac{\partial u\_i}{\partial t} + u\_j \frac{\partial u\_i}{\partial \mathbf{x}\_j} = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}\_j} + \nu \frac{\partial^2 u\_i}{\partial \mathbf{x}^2\_j},\tag{6}
$$

$$\nu\_t = \frac{\mathbb{C}\_{\mu}k^2}{\varepsilon},\tag{7}$$

where ν*<sup>t</sup>* is the turbulent viscosity, *C*<sup>μ</sup> is the constant. *k* and ε can be obtained from the following equations for turbulent kinetic energy and turbulent dissipation rate, respectively

$$\frac{\partial \mu\_j k}{\partial \mathbf{x}\_j} = \frac{\partial}{\partial \mathbf{x}\_j} \left\{ \left( \nu + \frac{\nu\_t}{\sigma\_k} \right) \frac{\partial k}{\partial \mathbf{x}\_j} \right\} + P\_k - \varepsilon\_\prime \tag{8}$$

$$\frac{\partial \boldsymbol{u}\_{\boldsymbol{j}} \boldsymbol{k}}{\partial \mathbf{x}\_{\boldsymbol{j}}} = \frac{\partial}{\partial \mathbf{x}\_{\boldsymbol{j}}} \left\{ (\mathbf{v} + \frac{\nu\_{\boldsymbol{l}}}{\sigma\_{\boldsymbol{c}}}) \frac{\partial \boldsymbol{\varepsilon}}{\partial \mathbf{x}\_{\boldsymbol{j}}} \right\} + \frac{\boldsymbol{\varepsilon}}{k\rho} (\mathbf{C}\_{\boldsymbol{c}1} \boldsymbol{P}\_{\boldsymbol{k}} - \mathbf{C}\_{\boldsymbol{c}2} \boldsymbol{\varepsilon})\_{\boldsymbol{\prime}} \tag{9}$$

Here, *C*ε1, *C*ε2, σε, and σ*<sup>k</sup>* are model constants. *Pk* is the turbulence generation term according to viscosity and buoyancy and is defined as follows:

$$P\_k = \nu\_t \frac{\partial u\_i}{\partial \mathbf{x}\_j} \left(\frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i}\right) \tag{10}$$

where the coefficients were defined as follows.

$$\mathbb{C}\_{\mu} = 0.09, \sigma\_k = 1.0, \mathbb{C}\_{\iota 1} = 1.44, \mathbb{C}\_{\iota 2} = 1.92,\tag{11}$$

### *5.2. Computational Domain and Grid System*

In this study, the total number of grid elements was 1,350,000. To investigate the change in the result value according to the number of grids, the grid elements were gradually increased from a small number. In the case of 551,000 grid systems, changes in flow distribution trends and internal velocity values in the fan and condenser were reduced (grid number 1,350,000), and changes in numerical calculation values, such as flow and velocity, became constant in more than 2,239,250 grid elements. Therefore, 1,350,000 grid systems were adopted in this study, as shown in Figure 5. Figure 6a,b shows the front and rear view of a computational grid system of the MCR. From Figure 6a,b, the arrangements of the compressor, condenser, and the CBM can be visualized. The internal velocity was measured between the compressor and condenser, as shown in Figure 6a.

**Figure 5.** Grid dependency of flow at MCR.

(**a**) Computational grid system (Front View) (**b**) Computational grid system (Rear View)

**Figure 6.** Computational grid system.

The governing equations were a type of continuously composed governing equations and were analyzed using the commercial CFD software ANSYS-CFX 16.2 [26,27]. The elements used in the calculation were a type of hexahedral elements, about 1,350,000, as shown in Figure 6. Moreover, a type of non-staggered grid system that determines pressure and velocity at different points was used to perform a linked calculation in pressure and velocity [5]. The CFX integrated coupled solver (pressure and velocity) was used as the solution method in this study.

### *5.3. Boundary Condition*

The discharge from the MCR was configured as an atmospheric condition, and the inflow by air through a fan was applied to the calculation by transforming the components measured using laser Doppler velocimetry (LDV), such as the velocity components of the fan along the axial, rotational, and radial directions of a cooling fan into X, Y, and Z components, respectively, in the Cartesian coordinate system [12,28]. The ideal air gas model was used in this study. The root mean square (RMS) convergence criterion was set to 10−4. The following equations are approximate equations for each test point using the velocity components calculated by using the LDV. The equation can be used to calculate the velocity for each point of the fan.

$$\text{V(a)} = -145232\text{X}^3 + 880.56\text{X}^2 + 11.54\text{X} - 0.35,\tag{12}$$

$$\mathbf{V}(\theta) = 45170\mathbf{X}^3 + 4474\mathbf{X}^2 - 38.64\mathbf{X} + 0.49,\tag{13}$$

$$\mathbf{V(r)} = \mathbf{64420X^3} - 1775.73X^2 - 29.55X + \text{ 0.69},\tag{14}$$

where V(a) is the velocity component along the axial direction, V(θ) is the velocity component along the circumferential direction, and V(r) is the velocity component along the radial direction, and they correspond to X, Y, and Z components, respectively, in the Cartesian coordinate system.

### **6. Experimental Results and Discussion**

### *6.1. Power Consumption and Inner Temperature of Chamber Room*

Figure 7 shows the power consumption for various temperature conditions of standard KS C 9305. From the Figure, it can be seen that the cause of an increase in power consumption was due to the increase in the chamber room temperature. As the temperature was raised by 1 ◦C, the power consumption was increased by 2 kW/month. The value of the power consumption can be calculated using the linear fitting function Equation (15).

$$\mathbf{y} = 2.06\mathbf{x} - 12.67,\tag{15}$$

**Figure 7.** Monthly power consumption for various temperatures in the chamber room.

### *6.2. Improvement of MCR Flow Line*

The characteristics of the MCR will cause heat transfer due to the forced convection if the flow line in the chamber is improved. By improving the flow line, it was possible to ensure more wind was received by the condenser by adjusting the flow line of the upper space of the condenser. In the improvement in the flow line, as shown in Figure 8, a box was attached to the upper wall of the MCR in which a 15 mm gap was left on the upper side of the condenser to avoid contact with it not to generate an insulation effect on its surface. Furthermore, the space between the compressor and the condenser where the flow finally passes through was configured as a curved surface to present a smooth flow.

(**a**) Original structure (**b**) Improved structure

**Figure 8.** Original and modified model of the machine room.

The temperature of the compressor and condenser installed inside the MCR for the original model was 53.9 ◦C and 43.2 ◦C, respectively. The temperature of the compressor and condenser for the improved model was 53.7 ◦C and 43.0 ◦C, respectively. Regarding the distribution of the average temperature of the specimens #1 and #2 inside the MCR in the improved model compared to that of the conventional model, the condenser showed a decrease in the temperature of 0.2 ◦C, and the compressor represented a decrease in the temperature of 0.2 ◦C. Thus, it represented an effect that decreased the monthly power consumption by 0.2 kWh, in which the reduction rate was 0.4%. By improving the internal flow line, it showed a decrease in the temperature of the whole surface of the condenser. Moreover, it is expected that the exit of the compressor that corresponds to discharge pipes also represented a decrease in the temperature. It is due to the fact that the flow, which had a high-speed element, progressed to the compressor via the condenser, except for the flow that intermittently passed through it in the resistance determined by the shape of the condenser by optimizing its upper space for the flow that progresses to the compressor. In addition, it was verified that the temperature in the condenser and compressor decreased due to the improvement of the flow line inside the MCR. Table 2 shows the reduction in power for the modified model when compared to that of the original model.

**Table 2.** Reduction in power of original and modified models.


### *6.3. Change in Position of the Cooling Fan*

The hot air in the MCR is discharged using a cooling fan. Thus, the airflow and internal flow rate can be varied according to the position of such a cooling fan, and this leads to present a difference in the heat transfer performance. As illustrated in Figure 9, the experiment was performed by adjusting the position of the cooling fan from a maximum of 40 mm to a minimum of 20 mm, with an interval of 5 mm referenced on the condenser.

**Figure 9.** Displacement of the cooling fan position.

In the experimental results of the monthly power consumption, as shown in Figure 10, the power consumption decreased when the position of the fan was closer to the condenser. Thus, it can be seen that the interval determined by 20 mm represented the highest efficiency. However, it is necessary to configure the position of the fan by considering the structure of the MCR as it is not possible to reduce the distance due to its structure excessively.

**Figure 10.** Monthly power consumption for various lengths between the condenser and the fan.

### *6.4. Optimization of Fan Flow Rate*

The airflow discharged from the fan inside the MCR becomes an essential factor for determining the performance of the whole heat transfer. In addition, the determination of an optimum condition for the fan operation is very important to improve the whole efficiency of the refrigerator because the power consumption operated by the fan according to the airflow represents large differences.

The airflow of the fan was obtained under the same conditions for a single product/system. The measurement method of the airflow for a single product/system was measured after adjusting motor revolutions per minute (RPM) of the motor. The RPM could be adjusted by controlling the frequency and voltage of the motor. This study applied a voltage-controlled adjusting method. The airflow of the fan was measured according to the position of each fan as an average value in

which 20 points were determined for measuring the airflow. Furthermore, the airflow converted as the measured value was presented, as shown in Figure 11.

**Figure 11.** Fan flow rate for various fan rotating speeds.

Figure 12 illustrates the results of the measurement of power consumption according to the airflow of the fan. The power consumption decreased as the airflow of the fan increased. It brought a decrease in the temperature of the compressor and condenser inside the MCR due to the increase in the airflow of the fan. Moreover, the decrease in such temperature increased the efficiency of the whole system. It can be seen that the results agree with the decrease in the temperature of the compressor and condenser as presented in the above two subsections that increased the efficiency of the refrigerator.

**Figure 12.** Power consumption for various flow rates of the cooling fan.

### **7. Numerical Results and Discussion**

### *7.1. Flow Patterns of Existing Refrigerator*

Figure 13a–c show the original CBM and the flow pattern of a model that applied conventional CBM. In the flow pattern, the internal air that cannot be discharged through the left vent of the compressor moved up to the upper section of the compressor and circulated to the upper section of the condenser. Then, the internal air was discharged to the upper left vent of the compressor according to the rise in the circulated air that moved up to the uppermost section.

(**c**) Velocity plots; top view

**Figure 13.** Original cover back machine (CBM) and flow pattern in the machine room of the refrigerator.

In the flow pattern at the space between the compressor and the condenser, the air passed from the condenser collided to the right side of the compressor, and that generated the circulated flow to the upper section of the condenser, and a whirlpool-like circulation was formed. As a result, it caused an increase in temperature at the entrance line of the condenser.

### *7.2. Selection of Modification Models*

In the results of the flow analysis for the existing models, it was necessary to ensure that the left vent for the flow that circulates at the left section of the compressor, the upper vent for the upper flow circulation, and the vent for the whirlpool and flow circulated due to the wall at the space between the compressor and the condenser. In addition, the vent was to be cut-off due to the internal flow from around the condenser to the compressor. The modified models were selected, as shown in Figure 14.

(**c**) Model 3

**Figure 14.** Modified structure of the CBM.

### *7.3. Numerical Simulation of the Modification Models*

Figure 15 illustrates the velocity field at a vertical section of +40 mm to the CBM referenced around the cooling fan. From the Figure, it can be seen that there was a smooth flow around the compressor in the modified model (c) due to the presence of a vent at the rear of the compressor when compared to models (a) and (b). In addition, it can be seen that the flow passed through the upper section of the compressor due to the fact that the flow stopped up at the vent of the condenser. In addition, model (d) shows lots of flows at the upper and lower sections of the compressor according to the installation of vents at the rear and around the compressor when compared with the model (c). In addition, there were smooth flows due to the discharge of flow that was generated at the space between the compressor and the condenser and was congested at the upper section of the compressor.

Figure 16 illustrates the velocity field at a vertical section of −20 mm to the CBM around the cooling fan. As shown in the figure, regarding the models (b–d) that ensure the vents at the left section of the compressor, it can be seen that the flow progressed around the compressor, whereas model (d) represents the most active movement. This shows that all modified models represented smooth flow around the compressor when compared to the original model.

**Figure 15.** Velocity plots of vertical section at +40 mm position from the center—front view.

**Figure 16.** Velocity plots of vertical section at −20 mm position from the center—front view.

Figure 17 illustrates the velocity field at a horizontal section of +20 mm to the CBM around the cooling fan. From the Figure, it can be seen that the models (b–d) showed the discharge of flows to the lower section of the compressor, whereas model (c) showed the flow that progressed to the compressor due to the cut-off of the vent of the condenser. Moreover, it can be seen that in model (d), flow was generated at the upper section due to the influence of the right upper vent of the compressor.

Figure 18 illustrates the velocity field at a horizontal section of −40 mm to the CBM around the cooling fan. From the Figure, it can be seen that the flow was discharged to the left section of the compressor. In the case of model (d), the flow was discharged as a vent was determined at the lower section of the compressor. In addition, because the horizontal sections in Figures 17 and 18 were not a normal symmetry but an axial symmetry, the sign of + to the horizontal center of the cooling fan and the—sign inside the MCR can be considered as the flow to the CBM.

**Figure 17.** Velocity plots of horizontal section at +20 mm position from the center—top view.

**Figure 18.** Velocity plots of horizontal section at −40 mm position from the center—top view.

Figure 19 and Table 3 represents a mass flow rate for the original and modified CBMs. Models 1 and 2 CBMs showed a high mass flow rate as the vent was generated at the left side of the compressor, and Model 3 CBM showed a high mass flow rate at the section of the compressor. Thus, it can be considered that the flow around the compressor is actively generated when compared to that of the existing models.

**Figure 19.** Mass flow rates of the original and modified CBMs.

**Table 3.** Mass flow rate (kg/s) of the original and modified cover back machine (CBM).


Since the modified models represented more smooth flows around the compressor when compared with the existing models, the decrease in the temperature at the compressor could be estimated. The temperature of the compressor and condenser installed inside the MCR for the original model was 54.5 ◦C and 43.5 ◦C, respectively. The temperature of the compressor and condenser for the improved model 3 were 52.0 ◦C and 42.7 ◦C, respectively. The temperature of the modified model 3 showed a decrease of about 2.5 ◦C and 0.8 ◦C at the compressor and condenser, respectively. As noted in Table 4, according to the decrease in the temperature, the power consumption of the refrigerator decreased by about 1.2% when compared to that of the existing models. Thus, the improvement in the flow line by modifying it in the MCR increases the efficiency of power consumption.


**Table 4.** Reduction in power for recent and modified models.

### **8. Conclusions**

An experimental and numerical study was carried out to analyze the performance improvement of a domestic refrigerator. In the results of the measurement of power consumption under the conditions of flow line control in the MCR and fan operation, the following conclusions can be summarized as follows:

The improvement of the flow line decreased the power consumption by 0.4% due to the decrease in the temperature of the compressor and condenser inside the MCR. The power consumption decreased in proportion to the distance between the fan inside the MCR and the condenser. Because the decreased power consumption was larger than the increase in the power consumed by the operation of the fan according to the airflow, the power consumption proportionally decreased according to the increase in the airflow.

For the modified Model 3 that improved the flow line around the compressor, it showed a smooth flow due to the lowest congested section. In the calculation of the mass flow rate for the existing and modified CBMs, the mass flow rate at the vent showed a high level in the model that vents around the compressor. In the smooth airflow inside the MCR in the modified CBM, it decreased power consumption by about 1.2% due to a decrease in the temperature of the section of the compressor (shell).

**Funding:** This research was funded by the Tongmyong University Research Grants 2017.

**Acknowledgments:** This research was supported by the Tongmyong University Research Grants 2017(2017A017). **Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2020 by the author. 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Investigation of Fluid-Structure Interaction Induced Bending for Elastic Flaps in a Cross Flow**

### **Tayyaba Bano 1, Franziska Hegner 2, Martin Heinrich <sup>1</sup> and Ruediger Schwarze 1,\***


Received: 30 July 2020; Accepted: 31 August 2020; Published: 5 September 2020

**Abstract:** With the recent increase in the design of light and flexible structures, numerical investigations of fluid and structure together play a significant role in most engineering applications. Therefore, the current study presents an examination of fluid-structure interaction involving flexible structures. The problem is numerically solved by a commercial software ANSYS-Workbench. Two-way coupled three-dimensional transient simulations are carried out for the flexible flaps of different thicknesses in glycerin for a laminar flow and Reynolds number ranging from 3 < Re < 12. The bending line of the flaps is compared with experimental data for different alignments of the flaps relative to the fluid flow. The study reports the computation of the maximum tip-deflection and deformation of flaps fixed at the bottom and mounted normal to the flow. Additionally, drag coefficients for flexible flaps are computed and flow regimes in the wake of the flaps are presented. As well, the study gives an understanding on how the fluid response changes as the structure deforms and the model is appropriate to predict the behavior of thick and comparatively thinner flaps. The results are sufficiently encouraging to consider the present model for analyzing turbulent flow processes against flexible objects.

**Keywords:** fluid-structure interaction; flexible flaps; bending line; deformation

### **1. Introduction**

Fluid-structure interaction (FSI) states the mutual interaction of a structure with a flowing fluid and develops the fluid and solid motions and the dependency of these motions on each other. Such kind of problems involve, the velocity and pressure of fluid effects and are affected by the structure deformation. In recent years, it has become a prominent characteristic of various engineering projects because design involves lighter and more flexible structures. The conduct of these structures with a fluid generates multi physical situations which are usually complex to handle. Numerical calculations play an important role in most engineering areas and have become more challenging due to the prediction of different flow regimes which cannot be seen in experiments (EXP) but it is still a challenge to have accurate numerical simulation of a FSI problem because of strong fluid-structure coupling. FSI problems are ever present in nature for example flapping wings, flags [1–9], aquatic animals [10], flexible vegetation [11–13], in hemodynamics [1,14] and human hair and hair sensors [15]. Thin structures are of great interest for microfluidic procedures in engineering and biomedical applications [16–18] and structural dynamics [18–20].

For investigating fluid motion together with the structure, extensive theoretical, experimental and numerical studies across flexible plates are carried out in the context of plate flapping, vibration, flutter instability, the effect of wake on plate stability and the estimation of natural frequency for

plates exposed to axial flow directions [21–25]. Results indicate that the flapping amplitude of the plate is mainly controlled by the elastic modulus and density ratio whereas frequency is a function of length-thickness ratio and density ratio. Also with an increase in Reynolds number, the plate indicates the transfer of symmetric to asymmetric flapping [21]. The forcing amplitude is the major influence on the plate resonance while flow velocity and plate bending rigidity indicate only the minor effect providing the forcing frequency is normalized with the natural frequency of the plate [10].

Aside from flat plates in axial flow direction, research work is offered to examine a vertical bottom-fixed flat plate in a laminar boundary layer. A study representing the Reynolds number, structure-to-fluid mass ratio and bending of the plate, respectively, determine the dynamic behavior of the flexible plate. Based on this, a two-dimensional (2D) numerical study [26] presents the dynamic response of a flexible plate engrossed in laminar condition for Reynolds numbers of 100, 200, 400 and 800 by an immersed boundary method (IBM). The dynamics of the plate indicate that the plate is stable for a low Reynolds number; however, when the Reynolds number reaches 800, the plate starts vibrating. The effect of bending rigidity, ranging from 0.1–1.6, demonstrates the periodic vibrations of the plate at low values and chaotic vibration when bending rigidity approaches between 0.2 and 0.4, then switching to periodic vibration again for a bending rigidity of 0.8 and 1.6. In comparison with the flapping flag [27], the plate shows a considerably different response in terms of frequency and vibrations.

The interaction of severely deformed structures with their environment is significant in many life science applications and technical systems. Nevertheless, precise and stable numerical methods for illustrating such problems are still exceptional. Tian et al. [28] developed a coupling between an existing immersed boundary flow solver together with a nonlinear finite element (FE) solver. Furthermore, they evaluate six practical applications of FSI specifically for three dimensions (3D), which primarily have large deformations. Considering a flexible plate as one of the test examples, the free-end deflection and drag coefficients measurements for Reynolds numbers ranging from 100–1600 with and without the effect of gravity and buoyancy forces are presented. Initially, the results of the baseline case with experimental data [29] are compared and then new simulations are offered for future benchmark cases. Moreover, the plate test case is used as a model to investigate the dynamic interaction of aquatic plants. The method overall leads to a numerical stability for all test cases.

The dynamic changes in the structure could possibly affect many aspects of fluid forces and present interesting fluid performance. Therefore, it is an important aspect for multiphase flow industrial applications. To compute such problems, tip-displacement of a rubber beam in a two-phase flow of air and oil was examined using a two-way coupled FSI [30]. A benchmark case [31] was simulated and previously published data of simulations and experiments were compared in terms of beam tip-deflection. The difference between 2D and 3D simulations has been clearly explained as well.

A scheme that has treated the deflection line of single and tandem beam configurations is presented by Axtmann et al. [32]. It is attained by interacting both experiments and direct numerical simulations (DNS) for Reynolds numbers ranging from 1–60 in a viscous fluid. Bending is calculated via beam bending theory. Bending lines of simulations and experiments demonstrate a good agreement for both single and tandem arrangements. A prediction model for cross flow velocity is offered. Moreover, local drag coefficient along pillar is computed and an empirical model based on drag coefficient is tested and proposed. Furthermore, a two-way coupled FSI study for different rows of flexible flaps is conducted to investigate the wavering behavior at the tip of flaps and development of vortices in the gap among the flaps. Experiments have been performed for water and glycerin against the silicon rubber flaps, reaching the maximum Reynolds number of 120, and the resulting tip-deflection of flaps is compared experimentally and numerically [5].

Investigating transient directional deformations of flexible obstacles by means of different numerical algorithms has been a subject of various studies involving the relevance of the extended arbitrary Langrangian Eulerian (ALE) method [33] and monolithic solvers [34,35]. The test problems presented [33] can be used as a benchmark for the formulation of numerical problems involving 2D FSI with immersed structures that have large displacements. While 2D benchmark cases, free-surface flows and test cases of an elastic obstacle and plate have been presented for monolithic approaches and validated with other numerical methods [34,35]. Similarly, for the analysis of hemodynamic problems, a new and more sophisticated computational formulation named immersed structural potential method (ISPM), based on the original IBM, is presented. The novel method accounts for the computation of FSI force field at the fluid mesh and description of cardiovascular tissue scheme by viscoelastic fiber-reinforced constitutive models. A new time-integration method is also introduced for the computation of a deformation gradient tensor. Three advancement cases of deformable cylinders and flapping membranes are presented. The study is applicable for highly complex FSI problems incorporating large deformations under pulsating flows and sinking and floating rigid bodies [36]. To describe the deflection of flaps and moving boundaries, the fluid- solid interface-tracking/interface-capturing technique (FSITICT) is applied to the FSI problem with ALE interface tracking and Eulerian interface tracking. The method is demonstrated by test FSI cases. First, test computations are performed for two flapping membranes and the second test case is composed of one flapping membrane and a movable, elastic arterial-structure. The flapping membranes and fluid are modelled in fully Eulerian coordinates in both test cases, while the elastic structure in the second case is computed in Lagrangian coordinates. The overall objective is to reduce the computational cost and test the applicability of the method for large deformations in order to analyze different cardiac cycles [37].

In the light of the above survey, it is concluded that the previous research efforts are made to encounter the dynamics of flapping structures in uniform flows. Prediction of the bending line and tip-displacement of a flexible flap fixed at the bottom in 3D is rarely seen in the literature. With the simple geometry of a flap in a rectangular fluid domain, the problem is fundamental in nature and can be used in applications that coincide with the flow, forcing functions and the solid boundary conditions. Therefore, this study emphasizes on the deformation characteristics of the flaps over time, flap bending positions and different flap orientations. The focus of this study is to explore the motion of an elastic flap immersed in a rectangular fluid domain and to analyze the bending of flexible structures in 3D under laminar flow conditions.

This paper is organized as follows-first, an insight into experimental procedure is presented then a computational model is demonstrated in terms of different flap thicknesses and positions. Later, the response of flexible flaps is analyzed in terms of bending and deformation in addition, the results of flap bending are offered in comparison with the experimental work. The effect of the different clamping angles is discussed as well. Similarly, the computation of drag coefficient at various clamping positions is a part of the current research. Flow regimes are revealed at various time steps and velocity contours in the wake are shown. Finally, some discussions are presented.

### **2. Material and Method**

A numerical two-way FSI model is presented for validating the towing-tank experimental data for flexible flaps. For the validation purposes, transient structure (flap) and unsteady flow are modeled using a commercial software ANSYS-Workbench. Each field is demonstrated separately by defining its geometry, properties, mesh, boundary conditions and respective solvers. Coupling is introduced between individual domains to establish a two-way FSI.

The experimental and numerical setups, settings and conditions are presented in the following subsets.

### *2.1. Experiments*

A summary of the experimental data is given in the present section, whereas detail can be found in the reference [32]. Flaps of different thicknesses are fixed in the center of a plate and immersed into the working fluid.

Flaps are designed with length of 100 mm, width of 20 mm and thicknesses of 5 and 10 mm. Several measurements for the deformable flaps are taken by varying the flow conditions (Reynolds

number, Re = 3 to 12) and flap arrangements that is, at clamping angles from 0◦ to 90◦ and inclined angle of 45◦. The bending of the flaps was recorded for each case. As the effective Reynolds number is small and viscous forces are dominating, only bending response is observed in the experiments instead of vortex shedding, which would cause periodic vibrations.

### *2.2. Numerical Model*

FSI modeling mainly deals with the information about how the fluid domain acts as a function of structure and how structural domain behaves as a function of fluid domain, at the interface, information is shared between the structural and fluid solver thus founding FSI [38]. In the present case, an elastic structure was surrounded by incompressible Newtonian fluid flow. Ωf(t) is defined as a time dependent fluid domain and Ωs(t) is defined as a time dependent structural domain. It is assumed that for all the time t, Ωf(t) is entirely occupied by the fluid and Ωs(t) is entirely occupied by the structure. Besides, Γo(t) = Ωf(t) ∩ Ωs(t) is an interface where the elastic structure interacts with the fluid. A schematic of FSI in terms of structural domain, fluid domain and interface for a present case of an elastic flap immersed in a rectangular fluid domain is shown in Figure 1.

**Figure 1.** A schematic of FSI in the present work.

For the established FSI solution, the kinematic, dynamic and geometric conditions at the common interface of fluid and structure should be satisfied, numerically at the interface Γo(t), the following conditions are assumed for all the time t:

$$\mathbf{v}\_i^s = \mathbf{v}\_i^f \tag{1}$$

$$
\sigma\_{ij}^s n\_i = \sigma\_{ij}^f n\_i. \tag{2}
$$

Here, *vs* and *v<sup>f</sup>* represent the structure and fluid velocity constituting the kinematic stability. Dynamic stability is attained by the continuous normal stresses of solid σ*<sup>s</sup>* and fluid σ*<sup>f</sup>* at the interface. *n* being the unit vector normal at the interface and the subscript *i* indicates its direction. According to the geometric condition stability fluid-solid domain always match and do not overlap [38].

Equations of motion for the structure and fluid are expressed as follows:

### **Structure equations**

The flexible flap is modeled as an elastic continuum therefore the balance of linear momentum is given by:

$$
\rho^s \stackrel{\circ}{\sigma^s\_i} = \sigma^s\_{i\bar{j},\bar{j}} \text{ in } \Omega\_\mathbb{s}(\mathfrak{t}).\tag{3}
$$

For linear elastic material, the linear Hook's law is followed by the structural stress as:

$$
\sigma\_{ij}^s = 2H\varepsilon\_{ij} + \lambda \varepsilon\_{kk} \delta\_{ij\prime} \tag{4}
$$

where ρ*<sup>s</sup>* is the structural density and σ*<sup>s</sup> ij*.*<sup>j</sup>* is the structural stress. The subscript *ij* of the stresses in Equations (3), (4), (6) and (7) represents the components of stress tensor whereas the subscript *j*, indicates the differential term involved in the linear elastic model.

*Appl. Sci.* **2020**, *10*, 6177

The structural stress is a function of strain ε and the Lamé parameters *H* and λ, which are defined as:

$$\begin{array}{c} \varepsilon\_{\bar{i}\bar{j}} = \frac{1}{2} (\mu\_{\bar{i},\bar{j}} + \mu\_{\bar{j},\bar{i}})\\ H = \frac{E}{2\left(1 + \mathbf{V}\right)}\\ \lambda = \frac{E\,\mathbf{V}}{\left(1 + \mathbf{V}\right)\left(1 - 2\,\mathbf{V}\right)} \end{array}$$

where *u* represents the displacement, *E* and Ю are Young's modulus and Poisson ratio, respectively [38–40].

### **Fluid equations**

The governing equations for time-dependent fluid regime are modeled by incompressible continuity and Navier-Stokes equations in terms of velocity *v<sup>f</sup>* , static pressure *P* and fluid density ρ*<sup>f</sup>* as follows [38–40]:

$$\operatorname{div} \dot{\boldsymbol{v}}\_{i}^{f} = \begin{array}{c} \text{0 in } \Omega\_{\mathbf{f}}(\mathbf{t}) \end{array} \tag{5}$$

$$
\rho^f \stackrel{\cdot}{\upsilon^f\_i} = \sigma^f\_{ij,j} \quad \text{in } \Omega\_\mathbf{f}(\mathbf{t}). \tag{6}
$$

The fluid stress for an incompressible fluid is simplified to:

$$\begin{aligned} \sigma\_{ij}^f &= -P\delta\_{ij} + \tau\_{ij} \\ \tau\_{ij} &= \frac{2}{3}\mu \Big(\varepsilon\_{ij} - \delta\_{ij}\mathfrak{c}\_{\mathbf{kk}}\Big) \\ \varepsilon\_{ij} &= \begin{pmatrix} \upsilon\_{i.j}^f + \upsilon\_{j.i}^f \end{pmatrix} .\end{aligned} \tag{7}$$

The dynamic mesh motion for an FSI problem is defined in terms of arbitrary movement of the grid points with respect to their frame of reference by considering the convection of the material points. The material derivative in terms of the movement of fluid mesh can be expressed as [40]:

$$\frac{d\upsilon^f}{dt} = \frac{\partial \upsilon^f}{\partial t} + \left(\upsilon^f - \mathsf{U}\right) \cdot \nabla \upsilon^f,\tag{8}$$

where *U* is the mesh movement of fluid, also *U* = *K Usur* representing *K* as a stiffness matrix of an elastic media and *Usur* being the fluid dynamics surface mesh movement.

Modeling of structure and fluid in ANSYS-Workbench are explained in the following sections:

### 2.2.1. Structural Model

The solid model is designed in 'Transient Structural' (ANSYS-Mechanical) by a (FE) analysis, afterwards the geometry is defined in the ANSYS Design Modeler. To estimate the flow as close as possible to the experiments, all dimensions and parameters are needed to be precisely defined according to the experimental setup and can be seen in Table 1.

Analogous to the experiments, the structural model is considered as a flexible flap that is, the bottom wall is fixed, while all other walls are free. The fixed wall is assigned as 'Fixed support' whereas the remaining part of the flap is given the 'Fluid Solid Interface.' The flap is meshed by hexahedral elements and are depicted in Figure 2a. Specifying 'w = 0.02 m' being the width of the flap and height as 5 w, the flap is positioned at 7.5 w (center) of the plate and 6.25 w from the side walls, whereas the distance of the flap from outlet is specified as 75 w. Accordingly, the rectangular flow domain is indicated as (150w, 20w, 12.5w). The structural model and FSI domain are shown in Figure 2a–c respectively. Part (a) of the figure represents the flap, its dimensions, mesh and transient conditions, (b) indicates the bending definitions and (c) shows the collected computational domain of fluid and structure.

**Figure 2.** (**a**) Flap configuration; (**b**) Bending definitions; (**c**) Computational domain.

### 2.2.2. Computational Model

The fluid domain is modeled by ANSYS-FLUENT for flaps immersed in a glycerin environment as depicted in Figure 2c. Because of the low Reynolds number range, a laminar fluid model is chosen. For an effective simulation in comparison to the experiments, a fixed coordinate system is used. The constant flow velocity as 'velocity inlet' is given for the input conditions, plate and flap are assigned as walls with a 'no slip' boundary condition while the top walls, remaining bottom parts and sides are assigned as 'slip walls' and 'pressure outlet' is defined for 0 Pa at the outlet conditions. A non-dimensional Reynolds number Re is defined as Re = *vw*/υ based on the flow velocity *v*, flap width *w* and fluid kinematic viscosity υ. The inlet flow velocities, *v*, ranging from 0.125–0.5 m/sec are given as input values based on the identified Reynolds number. The engineering parameters used for the computations are listed in Table 1. Due to subsonic incompressible flow field, pressure based solver with coupled scheme is used and second order implicit transient formulation is employed for modeling the unsteady flow.


**Table 1.** Engineering parameters and dimensions.

### **Five mm thick flap:**

Initially, simulations were performed for a flap of 5 mm thickness. To consider mesh independence, three different meshes were tested to explore the effect of mesh resolution on a flap bending. This mesh resolution was needed to compensate for the effect of computational time such as the computational time increasing after coupling the computational fluid dynamics (CFD) and FEA models [41]. Bending line convergence for different meshes was attained with coarser to finer: mesh 1: 52,486 elements, mesh 2: 101,641 elements and mesh 3: 137,693 elements. So, mesh 3 was chosen for the final solution, although this number of elements may vary since flap deforms and dynamic meshing takes place. Bending convergence is shown in Figure 3 whereas different meshes before and after flap deformations are depicted in Figure 4. Moreover, the difference between results obtained with mesh 2 and 3 is small, showing that grid convergence is attained. The difference between mesh 2 and 3 in terms of bending x and height y is 0.53% and 1.1% respectively.

**Figure 3.** Bending line comparison of different mesh resolutions at Re = 12.

As the large fluid domain is present, it is only desirable to have high mesh resolution around the flap and the dynamic mesh region in order to have influence of pressure and deformation gradient. For a 3D dynamic mesh domain in Fluent, a structured mesh cannot be applied [41]. For that reason, a tetrahedron unstructured mesh was chosen around the flap, though the rest of the fluid domain was designed as a structured mesh. To simulate the flap movement, both smoothing and remeshing methods were used in the dynamic mesh region such that the displacement of the flap was large compared to the grid size. Dynamic meshing is needed for both surface and volume mesh in 3D cases and it is necessary to retain element quality for both surface and volume meshes otherwise simulations may result in errors and computational differences [41]. Considering remeshing parameters, mesh scale information was followed and a minimum cell size of 6 <sup>×</sup> 10−<sup>4</sup> m and a maximum cell size of 0.04 m were used as input values. In addition, the interface is allocated as 'system coupling' in dynamic mesh zones under dynamic mesh settings with a cell thickness of 0.6 mm along the wall. To provide a balance

among stable mesh movement and total computational time, a time step Δt of 0.001 s is used for the coupled FSI problem for 2 s of physical flow time where the flap attained the steady-state condition.

**Figure 4.** (**a**) Different meshes used in the simulations, with increasing number of elements from left to right at initial time step; (**b**) deformed mesh 3 after 2 s.

### **Ten mm thick flap:**

In order to meet the experimental conditions, FSI analysis for a 10 mm thick flap is performed for different flap orientations, that is, varying the clamped angle from 0◦ to 90◦-for the angle of 0◦ the flap width is facing the flow, correspondingly angles are changed from 0◦–90◦. At 90◦, the flap is rotated completely and is parallel to the flow direction. Analysis is also performed for an inclined angle of 45◦. The 10 mm flap is modeled by the same structural mesh conditions as a 5 mm thick flap. Also for the fluid domain, the same mesh parameters are employed. Though the total number of elements vary in the range of 137,000–140,000 for different angles. The same boundary conditions and models are specified for 10 mm as described for the 5 mm thick flap. However, minimum and maximum length scales of remeshing differ, depending upon the flap positioning. Flap orientations designed in Design Modeler are shown in Figure 5a,b where (a) shows the different clamping positions for angle 0◦–90◦ and (b) indicates inclined position of the flap with the angle of 45◦.

**Figure 5.** (**a**) Flap positioning at different clamping angles angle 0◦–90◦; (**b**) Inclined angle.

### **3. Results and Discussions**

### *3.1. 5 mm Thick Flap*

A comparison between FSI and the experimental bending of a flap at different Reynolds numbers is expressed in Figure 6a–c. For maximum flow velocity of v = 0.5 m/sec (Re = 12) the flap bending is near to 1/2 of the flap length that is, y/L ≈ 0.5, for both FSI and experimental results. The simulated bending line represents the deflection of the flap in terms of x and y coordinates after reaching a steady-state. The flap bending increases with an increase in Reynolds number. However, a non-linear connection between flap deflection and Reynolds number is identified. Computational results can be taken for the whole flap length, though for the simplicity of measurements, five displacement points are measured along the length of the flap in experiments. In simulations, the x-coordinate bends to a maximum of x = 0.079 m and x = 0.076 m for leading edge and trailing edge respectively while in experiments, it is measured to be x = 0.070 m for leading edge and x = 0.067 m for the trailing edge. The results from experiments and simulations show the same bending behavior and are in good agreement. For the highest velocity, bending curves for both leading and trailing edge are shown in the Figure 6a,b, whereas the remaining results are focused only on leading edge of the flap.

Time-dependent total deformation curve of an elastic flap is processed by ANSYS-Transient Structural. Total deformation consists of directional deformations, d, in all three directions that is, X, Y and Z for a given flap. Resultant deformation varies with time from smallest to highest Reynolds number till it reaches its maximum value and then becomes constant. The deformation curves for all four Reynolds numbers at a physical flow time of 2 s after reaching steady-state can be seen in Figure 7 while Figure 8 indicates the flap visualization versus time till reaching steady-state for a Reynolds number of 9. The deformation is influenced by the magnitude of inflow velocity that is, maximum structural deformation occurs at the tip of the flap and has a least value for Reynolds 3 and increases by 0.0449 m reaching the peak value of 0.0982 m for a Reynolds number of 12. Besides, less

computational time is required by the flap for the higher velocity flow in converging to the steady-state value. Thus the FSI model consists of the observations for the maximum deflections offered by the flap.

**Figure 6.** Flap bending; (**a**) Re = 12 leading edge; (**b**) Re = 12 trailing edge; (**c**) Leading edge as a function of Reynolds number, thickness 5 mm.

**Figure 7.** Deformation verses time as a function of Reynolds number.

**Figure 8.** Flap visualization at time t, Re = 9.

### *3.2. 10 mm Thick Flap*

As described in Section 2.2.2, a flap of 10 mm thickness is used to simulate experimental data with different clamping orientations. For the highest Reynolds number of 12, the maximum deformation for a 10 mm thick flap is 0.53 times the 5 mm thickness for angle 0◦. However, it continues to decrease with increase in the clamping angle and shows the least at 90◦ which is 0.37 times the deformation at clamping angle of 0◦. Similarly, as the Reynolds number and the clamping angle increases, the x-coordinate of flap bending decreases indicating the minimum flap deflection with highest clamping angle (see Figure 9b) that is, the maximum x-coordinate for angle 0◦ is Δx = 0.039m whereas for 90◦ it bends to x = 0.011 m. Similarly, a small flap deflection is also observed for inclined angle of 45◦. Deformation at each position can be seen in Table 2 as well. Likewise, the comparison of experimental and computational results at various flap positions are shown in Figures 9 and 10 respectively.


**Table 2.** Deformation across flap orientations.

**Figure 9.** (**a**) Bending comparison at different thicknesses, Re = 9 clamping angle 0◦; (**b**) Bending verses Reynolds number for different clamping angles, thickness 10 mm.

**Figure 10.** *Cont.*

**Figure 10.** Inclined angle 45◦: (**a**) Flap bending; (**b**) Representation of flap deflection with respect to its initial state, thickness 10 mm

### *3.3. Drag Coe*ffi*cient*

For flow around solid bodies, the drag coefficient is defined as [42]:

$$\mathcal{C}\_D = \frac{F\_D}{\frac{1}{2}\rho v^2 S},$$

where *FD* is the drag force exerted by the fluid on the body and ρ is the fluid density, *S* is the body area perpendicular to the flow direction and *v* is the velocity of the flow. For the present simulations, the drag coefficient, *CD*, for three representative flexible flap cases is computed, that is, for a flap with a thickness of 10 mm and the clamping angles of 0◦, 90◦ and the inclined angle of 45◦ at a Reynolds number of 12. The frontal area of the flap facing the flow is defined for drag calculations meant for a current 3D case. Figure 11 represents the variation of drag coefficients versus time. The mean *CD* during t = 1.5–2 s for clamping angle of 90◦ is 3.61 and 1.95 for clamping position of 0◦. Whereas for an inclined angle of 45◦ it is smaller than the other two cases and gives a mean of 1.45. It appears that the more bended flap is a smaller obstacle to the flow and therefore reduces the resistance from the mainstream. Moreover, the drag coefficient values are comparable with a benchmark case of a flexible flat plate in a cross flow [28].

**Figure 11.** CD-comparison for different clamping positions.

### *3.4. Flow Regimes*

The velocity evaluation along with flap position for thickness 5 and 10 mm intended for clamped angle 0◦ and Reynolds 12 is displayed in the present section at various time steps. As discussed above (Section 3.2) and shown in Figure 12, the first two time steps at t = 0.05 s and 0.1 s corresponds to the initial deflection of the flap for both thicknesses. Nevertheless, the maximum deformation reached by the 10 mm thick flap is less for the 5 mm thick flap and attained in less time, whereas, for 5 mm thickness, the flap continues to bend and attain steady-state at approximately t = 1 s.

**Figure 12.** Flap movement of different thicknesses at various time steps.

The steady-state 3D configurations of clamping angle 0◦, 90◦ and inclined angle 45◦ at maximum inlet velocity is shown in Figure 13a–c presenting the region behind the flap where the transient bending converges to a steady-state and no flow fluctuation occurs downstream of the flap. The position of the flexible flap for the specified angles has changed from being straight to deformed. Flap wake indicates the absence of vortex shedding as observed in experiments because of a low Reynolds number range. Moreover, the flow past the flap is shown and compared for Reynolds numbers of 6 and 12 for clamping angle of 90◦ in Figure 13d,e and also for 5 mm thick flap representing the flow for minimum and maximum Reynolds number in Figure 13f,g, respectively. The higher Reynolds number leads to the higher velocity of the flow at the flap-tip and further downstream in the wake of the flap. Besides, the maximum bending of the flap in Figure 13g bends the streamlines near the flap downstream.

**Figure 13.** Comparison of flap wake; 10 mm Re = 12 (**a**) angle 0◦; (b) 90◦; (**c**) Inclined 45◦; (**d**) angle 90◦ Re = 6; (**e**) angle 90◦ Re = 12; (**f**) 5mm Re = 3; (**g**) Re = 12.

### **4. Discussion**

The present investigation focuses on the numerical investigation of flexible flaps in a glycerin environment for different thicknesses and flap orientations. The deformation of the flaps was measured to observe the various flap clamping positions under laminar flow and also the maximum tip-deflection of the flap was observed. For both experimental and simulation conditions, the flaps settle down to a steady-state and the bending patterns attained by current simulation model are very close to that observed in the experiments. The focus of the present FSI simulations is to design a 3D interaction between structure and fluid in order to narrow the literature gap for flexible structures. The results indicate that flap bending increases with the increase in Reynolds number in the range of

3 < Re < 12. However, for a different clamping position of 0◦ to 90◦ indicates minimum flap deflection with the highest clamping angle. Correspondingly, for an inclined position of 45◦ less deflection is perceived as compared to its undeformed shape. The objective is to judge the flap positions and flow regimes at varying computational time and to observe transient deformation response under different flow velocities. Likewise, the drag coefficients are also predicted for various flap alignments. Henceforward, the numerical model is suitable to accurately predict the response of thick and relatively thinner structures.

The flap bending in the current work is observed for low Reynolds number; however, the bending and periodic nature of flap motion is schedule by another 3D FSI model for higher Reynolds number in case of water as a working fluid. For this purpose, experiments will be performed for thinner flaps and the numerical model will be helpful in understanding the flow physics around the flaps. This will be further helpful in designing the flexible structures in hemodynamic, biomedical and structural dynamics.

**Author Contributions:** Conceptualization, T.B.; Methodology, T.B. and F.H.; Software, T.B.; Validation T.B.; Formal Analysis, T.B. and F.H.; Investigation, T.B.; Writing—original draft preparation, T.B.; writing—review and editing, T.B., F.H., M.H. und R.S.; supervision, M.H. und R.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded to T.B. by a scholarship from the funds of the Freistaat Sachsen Germany, grant number Landesstipendium and to F.H. by Deutsche Forschungsgemeinschaft (DFG) Germany, grant number DFG 1494/25-1.

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

### **Nomenclature**

The following symbols are used in the manuscript:


*Appl. Sci.* **2020**, *10*, 6177


### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

### **Numerical Simulation of the Oil Droplet Size Distribution Considering Coalescence and Breakup in Aero-Engine Bearing Chamber**

### **Fei Wang \*, Lin Wang \*, Guoding Chen and Donglei Zhu**

School of Mechanical Engineering, Northwestern Polytechnical University, Xi'an 710072, China; gdchen@nwpu.edu.cn (G.C.); zhudonglei@mail.nwpu.edu.cn (D.Z.) **\*** Correspondence: wfei@mail.nwpu.edu.cn (F.W.); wanglin@nwpu.edu.cn (L.W.); Tel.: +86-187-2900-1815 (F.W.)

Received: 12 June 2020; Accepted: 13 August 2020; Published: 14 August 2020

**Abstract:** In order to improve the inadequacy of the current research on oil droplet size distribution in aero-engine bearing chamber, the influence of oil droplet size distribution with the oil droplets coalescence and breakup is analyzed by using the computational fluid dynamics-population balance model (CFD-PBM). The Euler–Euler equation and population balance equation are solved in Fluent software. The distribution of the gas phase velocity field and the volume fraction of different oil droplet diameter at different time are obtained in the bearing chamber. Then, the influence of different initial oil droplet diameter, air, and oil mass flow on oil droplet size distribution is discussed. The result of numerical analysis is compared with the experiment in the literature to verify the feasibility and validity. The main results provide the following conclusions. At the initial stage, the coalescence of oil droplets plays a dominant role. Then, the breakup of larger diameter oil droplet appears. Finally, the oil droplet size distribution tends to be stable. The coalescence and breakup of oil droplet increases with the initial diameter of oil droplet and the air mass flow increasing, and the oil droplet size distribution changes significantly. With the oil mass flow increasing, the coalescence and breakup of oil droplet has little change and the variation of oil droplet size distribution is not obvious.

**Keywords:** bearing chamber; breakup; coalescence; population balance model; two phase flow; oil droplet size distribution; computational fluid dynamics

### **1. Introduction**

The bearing chamber is an important part in the aero-engine lubrication system. The lubricating oil is shed into the bearing chamber in the form of oil droplet via under-race lubrication method [1]. After the high-speed moving oil droplets collide with the wall, some of them are deposited to form the oil film, and others are splashed to form a large number of secondary oil droplets. Then, the complex air–oil two phase flow and heat transfer state of air, oil droplet, and oil film coexist in the bearing chamber. Because the oil droplet size is very small, many secondary oil droplets are suspended in the bearing chamber under the action of the air phase flow field. Furthermore, the interaction between oil droplets induces many dynamic events, such as collision, coalescence, deposition, evaporation, breakup, etc. The coalescence and breakup of oil droplets will affect the size distribution of secondary oil droplets and further restrict the research of flow and heat transfer in the bearing chamber. Moreover, the oil droplet size distribution directly affects the mixture state of oil and gas. Previously the Rosin–Rammler (R–R) distribution is widely used for the calculation of oil droplet size distribution in bearing chamber. However, experimental research indicates that R–R distribution is not consistent with the actual condition. Therefore, the study of droplet diameter distribution of the secondary oil droplets, which considers the coalescence and breakup of oil droplets, is of great significance in the aero-engine lubrication design.

At present, the research on coalescence and breakup of oil droplets mostly adopts the population balance model, which is earliest used in chemical process by Hulbert et al. [2]. Since then it has become a general method to analyze the particle diameter distribution of dispersed phase in multiphase system. The population balance model describes the time and space distribution of discrete phase by tracing the change of number density function of discrete phase with convection, diffusion and the source term of coalescence and breakup in the flow field. This method has been widely used because it remains part of the size information of discrete phase, and uses the averaging method to describe the motion of continuous phase, which reduces computation greatly. The core of population balance model lies in the profound understanding of coalescence and breakup mechanism of oil droplets. The coalescence and breakup of oil droplets is affected by the flow field around the droplets and the physical parameters (such as surface tension) of the oil droplet, which involves a very complicated multi-level and multiscale physical mechanism. The mechanism description of oil droplet breakup can be divided into three categories: (1) the droplet breakup caused by turbulence, (2) the droplet breakup caused by viscous shear force, and (3) the droplet breakup caused by the instability of the liquid–liquid interface. The coalescence mechanism of oil droplet is more complex. Moreover, the probability model is used more often, as the coalescence rate of oil droplet is equal to the product of collision frequency and coalescence efficiency between oil droplet. The collision mechanism of oil droplet can be divided into four categories: (1) collision caused by random turbulence fluctuation, (2) collision caused by buoyancy, (3) collision caused by shear force, and (4) collision induced by wave vortex. The coalescence efficiency model can be divided into three categories: energy model, velocity model, and liquid membrane drainage model.

During the past decades, many scholars have done many researches on the coalescence and breakup of bubbles/droplet and the size distribution of bubbles/droplet. Liao and Lucas et al. [3,4] reviewed the coalescence and breakup model of liquid, and introduced the research process of oil droplet coalescence and breakup mechanism. Ramkrishna et al. [5] investigated the effect of coalescence and breakup of oil droplet on the size distribution by using the population balance model, and further studied the flow behavior of liquid–liquid system from mechanism. Coulaloglou and Tavlarides et al. [6] assumed that the oscillating droplet will break up when the turbulent kinetic energy transmitted by the turbulent vortex to the droplet exceed the surface energy of the droplet, but the two droplets will coalesce when the contact time exceed the film discharge time of the two droplets. Tsouris and Tavlarides et al. [7] studied the droplet coalescence behavior in the turbulent using high-speed imaging technology and population balance model, proposed a novel collision frequency model, and compared the predicted droplet diameter with the experimental data. The results show that the model is applicable under different discrete phase fraction. Luo and Svendsen et al. [8] established a bubble breakup model, in which the bubble breakup rate is equal the product of the collision frequency with the turbulent vortex and the bubble breakup efficiency.

Regarding the solution of the population balance equation, Li et al. [9] used the quadrature method of moments (QMOM) to study the droplet breakup and coalescence behavior of the liquid–liquid dispersion system in the stirred reactor. The results show that the droplet coalescence behavior occurs with some difficultly in the low discrete phase fraction, due to the small probability of droplet collision. Li et al. [10] reviewed the CFD-PBM simulation in extraction column, and the principles and advantages and disadvantages of different solving methods of population balance model were compared. Xing et al. [11] used the experimental method and CFD-PBM coupled model to analyze the influence of viscosity in the bubble bed on the hydrodynamic behavior in the bubble area, such as the total gas holdup, the big and small bubble holdup, and the bubble size distribution. The results show that the CFD-PBM coupled model has a good prediction ability for the hydrodynamic behavior of the bubble bed in a wide range of gas velocities and viscosities. Wang et al. [12,13] used the population balance model to describe the bubble size distribution in the gas–liquid system, considered the different mechanism of bubble coalescence and breakup, established a relatively mature bubble coalescence and breakup model, and calculated the bubble size distribution in the gas–liquid system by numerical solution of

PBM. Zhao et al. [14] used a multi-Monte Carlo algorithm to simulate the collision and coalescence process of particles in nanoparticle flow, and the simulation results were in good agreement with the direct numerical simulation results. Chen et al. [15] analyzed the parameter relationship between the structural condition of bearing chamber and the oil droplet size distribution, clarified the oil droplet mass distribution based on the continuous oil drop size distribution, and introduced the concept of continuous oil droplet diameter into the analysis of oil droplet movement. Glahn et al. [16] measure the size distribution of oil droplet using the phase Doppler particle analyzer (PDPA) technology under various operating condition in the rotating disk chamber. In conclusion, the population balance theory has been widely applied in the analysis of coalescence and breakup between bubbles and droplets. However, no relevant research work has been done to consider the influence of coalescence and breakup on the oil droplet size distribution in the bearing chamber.

Therefore, based on the population balance model, the coalescence and breakup behaviors of oil droplets are described through the coalescence and breakup model of the oil droplet. Then, the air–oil two-phase flow model and the population balance model are solved by Fluent software. The numerical simulation of the coalescence and breakup processes of the oil droplet in a bearing chamber is performed. The air phase velocity distribution and the volume fraction of the oil droplets at different times are discussed. Moreover, the variation rules of oil droplet size distribution are obtained in different initial oil droplet diameter, air inflow, and oil inflow. The coalescence and breakup of oil droplets can change the size distribution of secondary oil droplets. Moreover, the oil droplet size distribution directly affects the mixture state of oil and gas. The research work in this paper can improve the accuracy of air–oil two-phase flow and heat transfer analysis, and provide reference for further research on oil mist concentration test in the aero-engine bearing chamber.

### **2. Numerical Model**

### *2.1. Geometric Model and Physical Parameters*

The geometric structure of the bearing chamber is shown in Figure 1. The structural parameters and values of the bearing chamber involved in the analysis mainly include the rotor radius of the bearing, *rs* = 60 mm; the height of the bearing chamber, *hb* = 30 mm; the width of the bearing chamber, *wb* = 30 mm; and the diameter of oil scavenge *ds* (scavenge is shown in Figure 2) and the air vent *dv* are both 16 mm. Seal air flows into the bearing chamber from air inlet and discharges through vent. The gap of oil inlet is 2 mm. Moreover, the lubricating oil is shed into the bearing chamber from the clearance between the bearing outer ring and the cage (oil inlet in Figure 1) and discharges through scavenge. The equivalent size of oil inlet is 4 mm (oil inlet in Figure 2).

**Figure 1.** The structural diagram of aero-engine bearing chamber.

In the analysis, the oil brand [17] is aviation lubricating oil 4109. The oil density ρ*<sup>l</sup>* is 926 kg/m3, the oil dynamic viscosity μ*<sup>l</sup>* is 0.007 Pa·s, and the oil surface tension coefficient σ*<sup>l</sup>* is 0.035 N/m. The gas density <sup>ρ</sup>*<sup>g</sup>* is 1.225 kg/m<sup>3</sup> and the air dynamic viscosity <sup>μ</sup>*<sup>g</sup>* is 1.789 <sup>×</sup> <sup>10</sup>−<sup>5</sup> Pa·s.

### *2.2. Mesh Model and Grid Independence Verification*

According to the structural parameters of the bearing chamber, Gambit software is used to build up the geometric model of the flow region, hexahedral grids are used to discrete the flow region, and the grid of air inlet and oil inlet is refined. The grid model is shown in Figure 2.

**Figure 2.** Grid model of the bearing chamber.

In order to eliminate the influence of mesh on the calculation, three groups mesh with different size are generated on the gambit platform, and the number of meshes is 487,585, 609,719, and 1,530,258, respectively. Under the condition of air inflow 12 g/s, the air phase velocities at four position are extracted on the section z = 0, and the four coordinate points are as follows; (−0.11,0,0), (−0.09,0,0), (0.09,0,0), (0.11,0,0). The calculation results are shown in Figure 3. From the figure, the velocities at four position obtained by three groups mesh are approximate. In order to save computing cost, the grid model of 609,719 is chosen in this paper.

**Figure 3.** Verification of grid independence, the number of meshes is 487,585, 609,719, and 1,530,258, respectively, when air inflow is 12 g/s, the air phase velocities at four positions (−0.11,0,0), (−0.09,0,0), (0.09,0,0), and (0.11,0,0) in z = 0 plane.

### *2.3. Mathematical Model*

In the air–oil two-phase flow, the dynamic events of oil droplet include collision, coalescence, breakup, deposition, etc. Under the combined effect of these dynamic events, the oil droplet diameter distribution changes with time. Under a series of physical actions, not only do the conservation of mass, momentum, and energy between oil droplet and air need to be considered, but also the population balance equation needs to describe the oil droplet diameter distribution. Especially, the characteristics of oil droplet collision, coalescence, and breakup are significant in the bearing chamber. It is necessary to focus on the change of oil diameter before and after coalescence and breakup. Therefore, the parameters, such as the volume fraction and velocity of each phase, are obtained by solving the Eulerian–Eulerian two-fluid model, and the dynamic events of oil droplet are described by using the population balance model. The oil droplet diameter distribution with the change of time and space is obtained by combining the two models. The Eulerian–Eulerian two-fluid model and the population balance model in this paper are as follows.

### 2.3.1. Two-Fluid Model

In two-fluid model [18], the air and oil phase are both treated as interconnected continuous media. As the volume occupied by one phase cannot be occupied by the other, the phase volume fraction is introduced. The volume fraction is a continuous function of space and time, and the sum of the volume fractions of each phase is 1. The volume fraction equation is as follows,

$$
\alpha\_{\emptyset} = \frac{V\_{\emptyset}}{V}, \; 0 \le \alpha\_{\emptyset} \le 1 \tag{1}
$$

$$\sum\_{q=1}^{n} \alpha\_q = 1 \tag{2}$$

where *V* is the total volume of air oil two-phase, *Vq* is the volume of the *q* phase, α*<sup>q</sup>* is the volume fraction of the *q* phase, and *n* has two phases: gas is labeled as phase 1 and oil droplet is labeled as phase 2. It should be noted that a two-fluid model can be considered a valid approximation only up to Stokes number ~0.1–0.2. A similar set of continuity equations and momentum equations can be derived for each phase.

*Appl. Sci.* **2020**, *10*, 5648

The continuity equation is

$$\frac{\partial(\alpha\_{\bar{q}}\rho\_{\bar{q}})}{\partial t} + \nabla \cdot (\alpha\_{\bar{q}}\rho\_{\bar{q}}\mathbf{U}\_{\bar{q}}) = 0 \tag{3}$$

where ρ*<sup>q</sup>* is the density of the *q* phase; **U***<sup>q</sup>* is the velocity vector of the *q* phase.

The momentum equation is

$$\begin{cases} \frac{\partial (a\_{\eta} \rho\_{\mathbf{i}} \mathbf{U}\_{\eta})}{\partial t} + \nabla \cdot (a\_{\eta} \rho\_{q} (\mathbf{U}\_{q} \otimes \mathbf{U}\_{q})) + \nabla \cdot (a\_{\eta} \mathbf{r}\_{q}) + \nabla \cdot (a\_{\eta} \mathbf{R}\_{q})\\ = -a\_{\eta} \nabla p + a\_{\eta} \rho\_{\eta} \mathbf{g} \pm \mathbf{M}\_{\eta} \end{cases} \tag{4}$$

where τ*<sup>q</sup>* is the viscous-stress tensors of the *q* phase, *p* is the pressure share by gas phase and the liquid phase, **R***q* is the Reynolds stress tensors, **g** is the gravitational acceleration vector, and **M***q* is the interfacial force term between the gas phase and liquid phase. When q is oil droplet, **M***<sup>q</sup>* is positive; when q is air, **M***<sup>q</sup>* is negative. The term **M***<sup>q</sup>* is usually composed by three forces: drag, virtual mass, and lift. In present work, only the drag force is taken into account, whereas the others are neglected. The drag force is calculated as follows,

$$\mathbf{M}\_{\emptyset} = \alpha\_1 \alpha\_2 (\frac{3}{4} \mathbb{C}\_D \frac{\rho\_l}{d\_{32}} |\mathbf{U}\_l|) \mathbf{U}\_r \tag{5}$$

where **U***r*= **U***<sup>g</sup>* − **U***<sup>l</sup>* is the slip velocity and **U***g*, **U***<sup>l</sup>* are the average vector of the air and oil, respectively. *CD* is the drag coefficient; *d32* is the mean Sauter diameter, calculated as the ratio between the moments of order three and two of the droplet size distribution.

In the two-fluid model, the turbulent Reynolds stress should be closed. The Reynolds stress arises from the operation of Reynolds averaging on the momentum equation. In this paper, the *k-*ε turbulence model [19,20] is selected to close the turbulence term in two fluid model. The turbulence model *k-*ε can be represented by the following equations,

$$\frac{\partial(a\_q k)}{\partial t} + \nabla \cdot (a\_q k \overline{u\_q}) = \nabla \cdot (a\_q \frac{\nu\_t}{\sigma\_k} \nabla k) + a\_q G\_k - a\_q \varepsilon \tag{6}$$

$$\frac{\partial(a\_{q}\varepsilon)}{\partial t} + \nabla \cdot (a\_{q}\varepsilon \overline{u\_{q}}) = \nabla \cdot (a\_{q}\frac{\nu\_{t}}{\sigma\_{\varepsilon}} \nabla \varepsilon) + a\_{q}\frac{\varepsilon}{k}(\mathbb{C}\_{k1}\mathbb{G}\_{k} - \mathbb{C}\_{k2}\varepsilon) \tag{7}$$

where *k* is turbulent kinetic energy; ε is turbulent energy dissipation term; ν*<sup>t</sup>* is turbulent viscosity, which can be expressed as *C*μ*k*2/ε, *Ck*1, *Ck*2, *C*μ, σ*k*, and σε are model constants; and *Gk* is the generating terms of turbulent kinetic energy.

### 2.3.2. Population Balance Model

The application of the population balance model [21] is used to calculate the oil droplet diameter distribution in the bearing chamber. The transport equation of the oil droplet number density function can be expressed as

$$\begin{aligned} \frac{\partial}{\partial t}[(n(V,t))] + \nabla \cdot \left[\overrightarrow{u}(n(V,t))\right] &= \underbrace{\frac{1}{2} \int\_{0}^{V} a(V - V', V') n(V - V', t) n(V', t) dV'}\_{1} \\ - \underbrace{\int\_{0}^{\infty} a(V, V') n(V, t) n(V', t) dV'}\_{\text{II}} &+ \underbrace{\int\_{0}^{\infty} b(V') \beta(V | V') n(V', t) dV'}\_{\text{III}} - \underbrace{b(V) n(V, t)}\_{\text{IV}} \end{aligned} \tag{8}$$

In the right-hand side of the equation in turn are the source term (I) generated by coalescence, the sink term (II) generated by coalescence, the source term (III) generated by breakup, and the sink term (IV) generated by breakup. Where *n*(*V*, *t*) is the number distribution function (1/m6), which means the number of oil droplets in the range of *V* and *V* + *dV* in unit volume is *n*(*V*, *t*)*dV*; *a*(*V*, *V*- ) is the coalescence rate function (m3/s) of oil droplets in volume *V* and *V* due to collision. *b*(*V*- ) is the breakup rate function (1/s) of oil droplets in volume *V*- . β(*V*|*V*- ) is the probability distribution density function (1/m3) of oil droplets in volume *V*, which breaks up from oil droplets in volume *V*- .

The population balance equation (PBE) is a hyperbolic integral-partial differential equation; only a few simple cases have an analytical solution, so it is necessary to use a numerical method to solve the PBE. At present, the usually used methods include the discrete method, moment method, and Monte Carlo method. Based on the calculation cost and accuracy, the discrete method is used to solve the population balance equation in this paper.

The basic idea of the discrete method is dividing the continuous diameter distribution of oil droplets into *N* discrete subintervals, and the size of all the droplets in the subinterval is equal to the node value *gi*. Meanwhile, the distribution function of the oil droplet number can be approximately as follows,

$$m(V, t) = \sum\_{i=1}^{N-1} N\_i(t)\delta(V - \mathfrak{g}\_k) \tag{9}$$

where *Ni* is the number of droplets in unit volume.

$$N\_i(t) = \int\_{V\_i}^{V\_{i+1}} n(V, t)dV\tag{10}$$

In most cases, the oil droplet size calculated by the coalescence and breakup model is not consistent with the node value of the subinterval, so it is necessary to allocate such droplets to the node value of the interval in a certain proportion, and the same time guarantee the conservation of the total oil droplet mass and the total number of oil droplet. Based on the above principles, Ramkishna proposes that the ratio of the newly generated droplets in the interval (*gi*, *gi*+1) to the node *gi* and *gi*+<sup>1</sup> is ψ(*V*, *gi*) and ζ(*V*, *gi*), and satisfies that

$$
\psi(V, \mathfrak{g}\_i)\mathfrak{g}\_i + \zeta(V, \mathfrak{g}\_{i+1})\mathfrak{g}\_{i+1} = V \tag{11}
$$

$$
\psi(V, \mathfrak{g}\_i) + \zeta(V, \mathfrak{g}\_{i+1}) = 1 \tag{12}
$$

By substituting Equations (9), (11), and (12) into Equation (8) and through a series of algebraic transformations, the discrete droplet transport equation is given by Equation (13).

$$\begin{aligned} \frac{d\mathbf{N}\_{i}(t)}{dt} + \nabla \cdot \left[ \overrightarrow{u^{\*}} N\_{i}(t) \right] &= \underbrace{\sum\_{j\_{r},k}^{|\mathcal{Z}| \ge k} \limits\_{j\_{r},k} \quad (1 - \frac{1}{2} \delta\_{j,\mathcal{k}}) x\_{i,j} a(\mathcal{g}\_{i}; \mathcal{g}\_{\mathcal{E}}) N\_{i}(t) N\_{k}(t)}\_{\mbox{\(\mathcal{G}\)} \leq (\mathcal{g}\_{i} + \mathcal{g}\_{\mathcal{E}}) \leq g\_{i+1}} \quad \text{(\(\mathcal{G}\)-1 and (\mathcal{G}\)-1)} \quad (13) \\ &\underbrace{-N\_{i}(t) \sum\_{k=1}^{M} a(\mathcal{g}\_{i}; \mathcal{g}\_{\mathcal{E}}) N\_{k}(t)}\_{\mbox{\(\mathcal{H}\)}} + \underbrace{\sum\_{i=1}^{M} \mathcal{Y}\_{i,k} b(\mathcal{g}\_{k}) N\_{k}(t)}\_{\mbox{\(\mathcal{H}\)}} \quad \text{(\(\mathcal{G}\)-1 and (\(\mathcal{G}\)-1)} \quad \text{(\(\mathcal{G}\)-1 and (\(\mathcal{G}\)-1))} \quad (14) \end{aligned}$$

where

$$\chi\_{i,j,k} = \begin{cases} \frac{\frac{\mathcal{S}\_{i+1} - V}{\mathcal{S}\_i + 1 - \mathcal{S}\_i}}{\frac{\mathcal{S}\_i - \mathcal{S}\_{i-1}}{\mathcal{S}\_i - \mathcal{S}\_{i-1}}} & V\_i \le V \le V\_{i-1} \\\ 0 & \text{otherwise} \end{cases} \tag{14}$$

After the droplet is broken up, the distribution proportion of each sub interval is expressed as

$$\gamma\_{i,k} = \int\_{\mathcal{S}^i}^{\mathcal{G}^{i+1}} \frac{\mathcal{G}^{i+1} - V}{\mathcal{g}^{i+1} - \mathcal{g}^i} \beta(V, \mathcal{g}\_k) dV + \int\_{\mathcal{G}^i}^{\mathcal{G}^{i+1}} \frac{V - \mathcal{g}\_{i-1}}{\mathcal{g}\_i - \mathcal{g}\_{i-1}} \beta(V, \mathcal{g}\_k) dV \tag{15}$$

Luo coalescence model [8] is adopted in this paper. The coalescence rate in model is defined as the frequency of new droplets generated after the collision of two droplets, and the volume of the two oil droplets is *Vi* and *Vj*. The coalescence rate is given by

$$\text{Ra}\begin{pmatrix} V\_{i\prime} \ V\_{j} \end{pmatrix} = \omega\_{\text{l\!\!\!\!\/\text\!}} \begin{pmatrix} V\_{i\prime} \ V\_{j} \end{pmatrix} \mathbf{P}\_{\text{l\!\!\!\!\/\text\!\!\/\text\!\!\/} V\_{i\prime} V\_{j} \end{pmatrix} \tag{16}$$

where <sup>ω</sup>*ag*(*Vi*, *Vi*) is the collision frequency; *Pag Vi*, *Vj* is the probability of coalescence. The collision frequency is obtained by

$$
\omega\_{\mathfrak{sl}}(V\_i, V\_j) = \frac{\pi}{4} (d\_i + d\_j)^2 n\_i n\_j \overline{u}\_{ij} \tag{17}
$$

where *uij* is the characteristic velocity of collision between two oil droplets with diameter *di* and *dj*, and the number density is *ni* and *nj*, which is defined as

$$
\overline{u}\_{ij} = \left(\overline{u}\_i^2 + \overline{u}\_j^2\right)^{0.5} \tag{18}
$$

$$
\overline{u}\_i = 1.43 \left( \varepsilon d\_j^2 \right)^{1/3} \tag{19}
$$

$$P\_{\rm af}(V\_{i\prime}, V\_i) = \exp\left\{-c\_i \frac{\left[0.75(1 + \mathbf{x}\_{ij}^2)(1 + \mathbf{x}\_{ij}^3)\right]^{0.5}}{\left(\rho\_2/\rho + 0.5\right)^{0.5}\left(1 + \mathbf{x}\_{ij}\right)^3} \mathcal{W}e\_{ij}^{1/2}\right\}\tag{20}$$

$$x\_{i\bar{j}} = \frac{d\_{\bar{i}}}{d\_{\bar{j}}} \tag{21}$$

$$\mathcal{W}c\_{i\bar{j}} = \frac{\rho\_i d\_i \left(\overline{u}\_{i\bar{j}}\right)^2}{\sigma} \tag{22}$$

where *ci* is constant, ρ*<sup>i</sup>* is the density of each phase.

Laakkonen breakup model [22] is adopted in this paper, which considers that the breakup of oil droplet is affected by both surface tension and viscous force; the degree of influence depends on the magnitude of two forces; and the collision model can be expressed as

$$b(V\_{\bar{j}}) = \mathbb{C}\_2 \varepsilon^{1/3} \text{erfc} \left( \sqrt{\mathbb{C}\_3 \frac{\sigma}{\rho\_l \varepsilon^{2/3} d\_{\bar{j}}^{5/3}} + \mathbb{C}\_4 \frac{\mu\_l}{\sqrt{\rho\_l \rho\_{\bar{\mathcal{S}}}} \varepsilon^{2/3} d\_{\bar{j}}^{5/3}}} \right) \tag{23}$$

where *C*<sup>2</sup> = 2.52, *C*<sup>3</sup> = 0.04, *C*<sup>4</sup> = 0.01.

The probability distribution function of sub droplets is given by

$$\beta \left( V\_i \middle| V\_j \right) = \frac{30}{V\_j} \left( \frac{V\_i}{V\_j} \right)^2 \left( 1 - \frac{V\_i}{V\_j} \right)^2 \tag{24}$$

The final breakup frequency is

$$
\Omega\_{\rm br} = b(V\_j)\beta(V\_i|V\_j) \tag{25}
$$

### *2.4. Initial and Boundary Conditions*

In the transient analysis, the boundary conditions include the inlet and outlet conditions of gas phase and lubricating oil, as well as the relevant wall conditions. The boundary conditions are treated as follows.


The calculation in this paper does not consider the effect of temperature. The oil droplets collision, coalescence, breakup, evaporation and deposition, nucleation and other dynamic events, as well as the geometry of the bearing chamber will affect the oil droplet size distribution in the bearing chamber. However, taking all the factors into account at the same time can make the analysis more complicated. Therefore, the calculation in this paper assumes that the temperature is constant and does not take into account the evaporation of oil droplets and the change of fluid properties.

### *2.5. Calculation Method*

The coalescence and breakup model of the oil droplet in the bearing chamber is solved by using ANSYS Fluent 17.0. The 3D unsteady model is adopted, the first-order upwind scheme is used, and the phase Coupled SIMPLE algorithm modified by pressure and velocity is employed for solving the transient equation. The coalescence and breakup model are as follows. The coalescence model uses the Luo model; the breakup model uses Laakkonen breakup model, in which both surface tension and viscous force are considered; and the oil droplet breakup is combine constrained by two forces, which is closer to the reality. The computation is performed on a Microsoft Windows 10 enterprise 64-bit operation system, and the simulations are run using 40 cores with 192 GB RAM machine.

### **3. Experimental Verification and Comparison**

In order to verify the correctness of the calculation method, the result in this paper is compared with the relevant test in literature [16]. In the literature, the Phase Doppler Particle Analyzer (PDPA) is used to sample 3000 oil droplet diameters in the aero-engine bearing chamber, and the histogram of oil diameter distribution is obtained. The distribution of droplet diameter is in a range of 14 to 120 μm. The operating condition in literature is air inflow m*<sup>g</sup>* = 10 g/s and oil inflow V*<sup>L</sup>* = 100 L/h. The oil droplet diameter distribution is calculated under these conditions. It should be noted here that in the calculation, the oil droplets larger than 120 μm in diameter are not considered, so only the oil droplet diameter distribution in a range of 10 to 120 μm is calculated. The calculation results are compared with the results in the literature in Figure 4. Through comparison, the calculation results are consistent with the experimental results in the literature, which verifies the rationality of the research method in this paper. The calculation results given in this paper can provide more accurate initial conditions for the analysis of the flow and heat transfer, and also can provide reference for the further test of oil mist concentration in the bearing chamber.

**Figure 4.** Comparison and verification of oil drop size distribution.

### **4. Results and Discussion**

In this paper, the distribution of air phase velocity field and the change of oil droplet volume fraction with time in the bearing chamber at different times are analyzed when the air inflow is 12 g/s, oil inflow is 15 g/s and initial oil droplet diameter is 28 μm. Moreover, the influences of different initial oil droplet diameter, air inflow, and oil inflow on the oil droplet diameter distribution are further analyzed. The operating condition parameters are shown in Table 1.



According to the mechanism of oil droplet coalescence and breakup, the turbulence of air phase flow field has a great influence on the coalescence and breakup of oil droplet. Therefore, the velocity of air phase in the bearing chamber is firstly analyzed. Figure 5 shows the distribution of air phase velocity in the bearing chamber at different time when the air inflow is 12 g/s, oil inflow is 15 g/s, and the initial oil droplet diameter is 28 μm. As can be seen from the figure, the velocity near the outlet of the bearing chamber is largest, and the other region is smaller. Because the seal air enters into the bearing chamber from the air inlet, it first moves along the axial direction. When the air is blocked by the wall surface, the velocity changes, which results in an uneven velocity distribution of the flow field. The air velocity in the middle region is significantly lower than that in the inner wall surface and the outer wall surface. Meanwhile, many vortices with different size existed in the bearing chamber. With the passage of time, the flow field structure tends to stable.

**Figure 5.** The velocity (m/s) distribution of air at different time when the air inflow is 12 g/s, oil inflow is 15 g/s, and the initial oil droplet diameter is 28 μm at z = 0 plane (the far left is vent; the far right is scavenge). (**a**) *t* = 50 ms; (**b**) *t* = 100 ms; (**c**) *t* = 200 ms; (**d**) *t* = 400 ms.

Figures 6 and 7 show the volume fraction of oil diameter in the range of 56 to 80 μm and 28 to 40 μm when the air inflow is 12 g/s, oil inflow is 15 g/s, and initial oil droplet diameter is 28 μm, respectively. In Figure 6, the volume fraction of the large diameter oil droplet on the outer wall surface of the bearing chamber is large, that near the inlet is relatively low, and gradually increases along the z axial direction. The volume fraction of large diameter oil droplet near the outer wall of the bearing chamber increases gradually with the passage of time. The reason is that the oil droplet flow with the air phase, and enter the bearing chamber from inlet. During the flow process, the oil droplets with smaller diameter constantly gather. When the air phase flow field reaches the wall that is opposite to the inlet, the flow direction of air phase changes and flows along the wall. During the movement of the oil droplet with the air phase flow field, the oil droplets constantly gather again. Therefore, the volume fraction of oil droplets at the external wall surface of the bearing chamber is large. It can also be seen from the figure that when the oil droplet with larger diameter runs to the entrance with larger velocity, the breakup will occur. The oil droplets are coalesced and broken in the bearing chamber, and finally reach a steady state. In Figure 7, the small diameter oil droplet is mainly distributed near the inlet, and this distribution gradually decreases along the axial direction. It can be seen that the size of the oil droplet in the bearing chamber is different because of the coalescence of oil droplets. The small diameter oil droplet is mainly concentrated near the inlet, while the large diameter oil droplet is mainly concentrated on the outer wall of the bearing chamber.

**Figure 6.** Oil droplet volume fraction (%) at different time with diameter in the range of 56 to 80 μm when the air inflow is 12 g/s, oil inflow is 15 g/s, and initial oil droplet diameter is 28 μm (the oil inlet is on the bottom, the vent is in front and the z axial direction is up). (**a**) *t* = 50 ms; (**b**) *t* = 100 ms; (**c**) *t* = 200 ms; (**d**) *t* = 400 ms.

**Figure 7.** Oil droplet volume fraction (%) at different time with diameter in the range of 28 to 40 μm when the air inflow is 12 g/s, oil inflow is 15 g/s, and initial oil droplet diameter is 28 μm (the oil inlet is on the bottom, the vent is in front and the z axial direction is up). (**a**) *t* = 50 ms; (**b**) *t* = 100 ms; (**c**) *t* = 200 ms; (**d**) *t* = 400 ms.

XUH

### *4.2. The Change Rule of Droplet Size Distribution with Time*

The change of oil droplet diameter distribution in the bearing chamber with time when the air inflow is 12 g/s, the oil inflow is 15 g/s, and the initial oil droplet diameter is 28 μm is shown in Figure 8. The oil droplets enter into the bearing chamber, due to the flow field development is not sufficient oil droplets coalesce firstly in the initial stage. Moreover, the volume fraction of large diameter droplets increases simultaneously. Then, the breakup of the oil droplet emerges prominently under the action of gas phase turbulence, and the volume fraction of big diameter droplets decreases. At last oil droplet diameter distribution tends to stable and the size of oil droplet mainly concentrates in the vicinity of 50 μm and 80 μm. From another angle, it shows that the oil droplet size distribution in the bearing chamber is not a simple R–R distribution

**Figure 8.** Droplet size distribution at different time when the air inflow is 12 g/s, the oil inflow is 15 g/s, and the initial oil droplet diameter is 28 μm.

### *4.3. The Influence of Di*ff*erence Initial Diameter of Oil Droplet on the Size Distribution*

Figure 9 shows the influence of different initial diameter on the oil droplet size distribution in the bearing chamber when the air inflow is 12 g/s and the oil inflow is 15 g/s at *t* = 250 ms. It can be seen that the influence of the initial diameter of the oil droplet on the volume fraction of different diameter oil droplet is obvious, and the peak diameter moves in the direction of large diameter. When the initial diameter increases to 56 μm, the peak value of oil droplet diameter increases to 70–80 μm. The size of the oil droplet diameter determines the value of the inertia force. The small diameter oil droplet is greatly affected by the air flow field, while the large diameter oil droplet is the opposite. The increasing of the inertia force of oil droplet weakens the moving velocity of the oil droplet, so the coalescence of oil droplet strengthens. Therefore, the proportion of large oil droplet is larger with the initial diameter of oil droplet increasing in the bearing chamber.

**Figure 9.** The influence of initial diameter on oil droplet size distribution when the air inflow is 12 g/s and the oil inflow is 15 g/s at *t* = 250 ms.

### *4.4. The Influence of Air Inflow on the Oil Droplet Size Distribution*

The change of oil droplet size distribution under different air inflow is shown in Figure 10; when the oil inflow is 15 g/s, the initial oil droplet diameter is 28 μm at *t* = 250 ms. It can be seen that the peak value of oil droplet diameter is mainly concentrated around 50–80 μm. When the air flow is small, the proportion of oil droplet with 50 μm is larger than that with diameter of 80 μm. With air inflow increasing, the effect of turbulence intensity increases, and the strong disturbance in the flow field and the entrainment of vortices enhance the mutual collision between oil droplets. After the oil droplets collide with each other, the proportion of large diameter oil droplets increases in the oil droplet group due to the adhesion. Meanwhile, the increase of air inflow leads to the increase of oil droplet velocity, which promotes the coalescence between oil droplets. The smaller diameter oil droplet coalesces into larger diameter oil droplets. When the air inflow is 20 g/s, the peak value of oil droplet diameter distribution is mainly around 80 μm. This shows that the air inflow has a significant effect on the oil droplet diameter distribution.

**Figure 10.** The influence of air inflow on the oil droplet size distribution when the oil inflow is 15 g/s, the initial oil droplet diameter is 28 μm at *t* = 250 ms.

### *4.5. The Influence of Oil Inflow on the Oil Droplet Size Distribution*

The oil droplet diameter distribution under different oil inflow when the air inflow is 12 g/s, the initial oil droplet diameter is 28 μm at *t* = 250 ms is shown in Figure 11. It can be seen that the peak value of oil droplet diameter distribution is concentrated at 45–80 μm. The effect of oil inflow on the distribution is not obvious because the increasing oil inflow only increases the volume fraction of oil in the bearing chamber, but does not change the size distribution of the oil droplets.

**Figure 11.** The influence of oil inflow on the oil droplet diameter distribution when the air inflow is 12 g/s, the initial oil droplet diameter is 28 μm at *t* = 250 ms.

### **5. Conclusions**


The oil droplet size distribution directly affects the mixture state of oil and gas in the bearing chamber. Compared with the R–R distribution, the oil droplet size distribution obtained in this paper is more consistent with the real situation by considering the influence of coalescence and breakup of oil droplets. It can provide more accurate initial conditions for the future calculation of heat transfer coefficient internal and external the bearing chamber.

**Author Contributions:** Conceptualization, F.W.; Methodology, F.W.; Software, F.W.; Validation, F.W.; Formal Analysis, F.W. Investigation, F.W. and L.W. Writing—original draft preparation, F.W.; Writing—review and editing, F.W., L.W., G.C., and D.Z.; Supervision, G.C.; Project Administration, G.C.; Funding Acquisition, L.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 51975475, and the Fundamental Research Funds for the Central Universities, grant number 31020200503002.

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

### **Nomenclature**


### **Abbreviations**


### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

*Article*

### **A New Vector-Based Signal Processing Method of Four-Sensor Probe for Measuring Local Gas–Liquid Two-Phase Flow Parameters Together with Its Assessment against One Bubbly Flow**

### **Xiaohang Qu 1, Qianjian Guo 2,\*, Yi Zhang 1,3,\*, Xiaoni Qi <sup>1</sup> and Lei Liu <sup>3</sup>**


Received: 15 July 2020; Accepted: 5 August 2020; Published: 7 August 2020

**Abstract:** A multiphase flow measurement technique plays a critical role in the studies of heat and mass transfer characteristics and mechanism of the gas–liquid two-phase, the practical measurement of the gas–liquid flow and the improvement of multiphase theoretical models. The four-sensor electrical probe as an emerging measurement method has been proved to be able to get the local flow parameters of multi-dimensional two-phase flow. However, few studies have been reported using the four-sensor probe to obtain the interface information (e.g., the interface direction and velocity). This paper presents a new signal processing method by which the interface direction and velocity can be obtained, besides void fraction, interfacial area concentration (IAC) and bubble chord length. The key solution is to employ the vector-based calculating method, which possesses the merits of simplicity and efficiency, to gain the interface velocity vector through legitimately assuming a direction of the interface velocity. A miniaturized four-sensor electrical probe was made and a gas–liquid two-phase flow experiment was performed to test the proposed signal process scheme. The two-phase flow was controlled to be in cap-bubble flow regime. To validate the availability and reliability of the proposed method, the local flow parameters obtained by the probe measurement were compared with the results from visual measurement technique in the same flow conditions. The comparison indicates that the above local flow parameters from four-sensor probe measurement are in good agreement with the visual measurement results, with maximum deviations of chord length of 8.7%, thereby proving the correctness of the proposed method.

**Keywords:** gas–liquid two-phase flow; local flow parameters measurement; four-sensor probe; vector-based signal processing; visual measurement experiment

### **1. Introduction**

Gas–liquid two-phase flow is a common phenomenon occurring in petroleum, chemical, refrigeration and power generation industries [1–3]. Due to the unstable flow, heat and mass transfer process, the flow pattern and interface structure of two-phase flow are usually complex [4–6]. Therefore, the two-phase flow parameters are difficult to measure quickly and accurately [7]. However, the knowledge of local parameters of two-phase flow plays a critical role in the studies of the heat and mass transfer characteristics and mechanism of the gas–liquid two-phase, the development of two-phasemodeling research, the optimization of two-phase flow patterns and the safety and stability of equipment operation [8–10]. The theoretical models, including the two-fluid model and interfacial area transportation model, rely heavily on the

advanced measurement techniques to provide benchmark database and possibility of inspecting new phenomenon and physical laws. The flow parameters including phase distribution, void fraction, interfacial area concentration (IAC), bubble size, velocity, etc. of two-phase flow have great impacts on the heat and mass transfer characteristics, the reaction efficiency and operation safety of multiple chemical applications [11–15]. Hence, it is particularly necessary to develop and utilize accurate, fast and convenient methods to explore the two-phase flow details.

Currently, the measurement methods of local two-phase flow parameters can be divided into the following two categories: (1) the photography and image process techniques, having the advantage of lack influence on the flow, usually called visual measurement [16–18] and (2) the point by point measurement [3,19,20], in which optical or electrical signal that can be altered by involved phases. However, the first kind of methods is only applied to systems where the flow channel or vessel is transparent, unless the photography employs high energy ray to penetrate metal walls [21,22]. Besides, it is hard for visual measurements to distinguish bubbles at different depths when bubbles are overlapped in complex two-phase flow [23]. Nevertheless, this method has contributed greatly in the identification of flow regimes, track of simple bubbly flow and determination of bubble condensation rate in phase changing flows [23–28]. For the second category of methods, one of its important kinds is the wire-sensor mesh [29–33] where a matrix of measurement points is created at the cross points of two arrays of parallel mesh wires. Electrical signals from the wire-sensor mesh are collected and then analyzed by a program to get void fraction, bubble size, bubble velocity and so on across the whole flow field. Another kind is the sensor probe including single-sensor probe, double-sensor probe, four-sensor probe and other multiple-sensor probes [34–42]. Moving the probe across the two-phase flow field, local flow parameters at different positions can be measured. Although both the wire-sensor mesh and sensor probe interfere measurements affect the original two-phase flow field, this intrusion can be limited by a great extent with very thin wires or probes. Moreover, compared with the visual measurement, the point-to-point sensors are easier to be adopted online for practical opaque tube and vessel and more complex two-phase flow.

In the above-mentioned method, only two sensors apart from the probe shell are used in the double-sensor probe, and thus two series of signals can be collected. In the simple measurement theorem of double-sensor probe, it is assumed that the direction of interface and velocity of bubbles (or interface) are both in parallel to the connection between the two sensors' tips. As a result, the double-sensor probe can be applied to measure local flow void fraction, IAC and magnitude of velocity when the two-phase flow is very consistently stable bubbly flow, and the accuracy can be very high without too much care to the signal process and correction algorithm [34]. However, more assumptions and more complicated signal processing and correction algorithms are necessary if the flow becomes complicated or more parameters, e.g., bubble size and interface direction, have to be detected. Inevitably, more uncertainties would be produced in the application of single-sensor probe or double-sensor probe if employing these assumptions [43–46]. As uncertain assumptions must be incorporated, it is hard for the double-sensor probe to obtain real local flow parameters in multi-dimensional two-phase flow.

Meanwhile, the typical multiple-sensor probe, namely the four-sensor probe, as a promising alternative compared with single-sensor probe and double-sensor probe, is developed by a few researchers to omit the involvement of so many assumptions. According to measurement theorem of four-sensor probe, it is unnecessary to coordinate the connections of the sensors' tips to the direction of interface or direction of interface displacement velocity, while the void fraction, IAC and chord length can still be calculated. In addition, the interface direction can be obtained by the proper processing and trigonometric operation of the four series of signals from the four sensors [47–51]. However, the bubble velocity or interface displacement velocity is still unknown, except its component in the normal direction of interface [36]. It has been proved that this problem can be compensated by employing several four-sensor probes or multiple-sensor probe, five-sensor probe, six-sensor probe, etc. [52,53]. With a much simpler assumption of the two-phase flow, obtaining the full velocity vector by using merely one four-sensor probe can also be realized. For instance, assuming a spherical or symmetric

bubble shape and a bubble velocity perpendicular to the symmetric plane has been demonstrated to be a concise way to get the full components of bubble velocity [54–56]. The above simplifications when calculating the total interface velocity vector only apply to special cases. Considering this fact, the local interface velocity directions can be practically provided by prior measurement, flow simulation and even other legitimate assumption in many two-phase flows. It is necessary and significant for obtaining interface details to propose a new method to get the total interface velocity vector based on the known interface velocity directions.

Besides the development of algorithm for four-sensor probe, the miniaturization of the probe and sensors are essential to guarantee the accuracy of measurement. According to the published literature [47,50], both the diameter of sensor wire and the total front area formed by all four sensors have great influence on the bubbles' behavior. Thus, it is recommended to use as small as possible sensors and probe of cross-sectional area. With the assumption of a much smaller probe than the size of bubble, the error produced during the application of the probe stems mainly from the bouncing away of bubble from the probe and the slipping away of interface through the gap between sensors. Their errors or uncertainties are acceptable [36,47,49]. Therefore, with the solution of the many issues of four-sensor probe, including its fabrication and further miniaturization, the correction for its disturbance to original flow and the improvement of signal process algorithm, it is promising to be widely applied in multiphase flow measurement inside chemical reactor, oil piping, power generation facility and heat exchanger. In addition, the method will be easy to be adopted for different combinations of fluid components. Even if the flow were experiencing phase changing, including boiling and condensation, the method would still be applicable with proper algorithm improvement and correction.

However, the bubbles in two-phase flow keep deforming and are hard to be taken as sphere or symmetric, resulting in few studies on the accurate, convenient and efficient measurement of local flow parameters having been published thus far. Therefore, a new signal processing method for four-sensor probe to get the bubble velocity vector was developed in this study and the vector-based calculation was used for the first time to deduct the local flow parameters. Besides, the interface direction obtained from the probe were for the first time validated against a visual experiment that was also performed. The void fraction, IAC, bubble velocity, bubble chord length and interface direction resulting from the probe measurement were compared with the visual measurement. The application perspectives of this method in the field of mass and heat transfer of gas–liquid two-phase flow is discussed. The proposed methods in the present paper are expected to be useful in the heat and mass transfer characteristics and mechanism studies of the gas–liquid two-phase and the direct measurement of two-phase flow. Meanwhile, it can able provide significant database for the improvement of two-phase models.

### **2. Measurement Principles of Four-Sensor Probe**

### *2.1. Electrical Circuit of Four-Sensor Probe*

As shown in Figure 1, an electrical circuit can be adopted in the four-sensor probe measurement. Since signal filtering and noise reduction can be easily realized by the appropriate signal pre-processing MATLAB codes (MATLAB 2017b, MathWorks, Inc., Natic, MA, USA), the circuit elements responsible for these functions were not necessary and thereby only four electrical resistances were employed in the circuit. The probe contains four sensors denoted *s*0, *s*1, *s*<sup>2</sup> and *s*3, which connect to the negative electrode of the DC power supply through their respective resistance. The four sensors are covered by a rigid stainless-steel shell, which connects to the positive electrode of the DC power and functions as the common high voltage pole.

The shell of the probe made of metal was not insulated and thus always in contact with water. The four sensors are well insulated except at their very ends where the sensors can contact the water and thus get through to the positive pole of DC supplier. The signals from *s*<sup>0</sup> to *s*<sup>3</sup> are either high level or low level, depending on whether the sensor tips are submerged into water or exposed to the bubble air. A high voltage level indicates that a particular sensor is submerged in water and a low voltage

level indicates that a particular sensor is exposed to air. Finally, the signals of high or low voltage from four sensors were transmitted to data collection system composed by a data acquisition unit (Art Technology, Beijing, China) and a PC. It is worth noting that the sensors have to be connected to the negative electrode of the DC power to avoid electrochemical corrosion and expand the probe's lifespan.

**Figure 1.** Schematic diagram for the measurement circuit of the four-sensor probe.

For the data acquisition system, it is suggested in the literature that a sample frequency higher than 10 kHz is required to guarantee resolution. Therefore, a sample frequency of 10 kHz for each of the four sensors was employed in this study.

### *2.2. Fabrication of the Probe*

As the impact of probe on bubble only becomes negligible with small probe size, the accuracy of the measurement improves with the downsizing of the four-sensor probe, thus the size of the probe should be miniaturized as much as possible. The four-sensor probe used in the present study was hand-made using stainless steel tube with internal diameter 1.2 mm and outer diameter 1.5 mm as the positive electrode and four copper wires of diameter 0.1 mm as the negative electrodes. The specific dimensions of the probe and its picture in reality are shown in Figure 2.

**Figure 2.** Dimensions of the four-sensor probe. (**a**) Probe in side view; (**b**) Probe in front view; (**c**) Probe photo.

The dimensions of the four-sensor probe and the relative position of sensor tips are shown in Figure 2a,b, respectively. As shown in Figure 2c, to provide support for the probe, the four sensors made of copper wire were fixed inside a stainless tube by epoxy resin, and the stainless tube with a length of 300 mm was then fixed by its end inside a short tube of internal dimeter of 5 mm with resin, before mounting it in a slide module, as shown in Section 3.1. Each copper wire was covered by electrically insulation material except the tip was polished by sandpaper for electrical contact with water. Choosing the wire tip as coordinate center and the axis of *s*<sup>0</sup> as the *y*-coordinate, the three vectors formed by the probe are *S*<sup>1</sup> (0, 1.75, 0.5), *S*<sup>2</sup> (0.5, 1.5, 0.5) and *S*<sup>3</sup> (0.5, 1.5, 0), which are used in Section 2.4.

### *2.3. Signal Pre-Processing*

After data collection, the signals were processed by a MATLAB program consisting of mainly two functions, the pre-processing and deduction of local flow parameters. The signal pre-processing was designed to obtain the time instants when each of the four sensors penetrate or recede a bubble.

In an ideal case, the signals should be square wave (see Figure 3d) with the high-level representing sensor contact with water and the low-level representing passing-by of a particular bubble. However, mainly due to the delaying of data collection system and electromagnetic interference, the practically collected signal demonstrates noising inclining and fluctuating features, as shown in Figure 3a. Average filtering was applied to attenuate the noise signal firstly, with the result shown in Figure 3b. Then, a threshold voltage was chosen for signal binarization, as shown in Figure 3c. The signal inversion so that the high level corresponds to air is shown in Figure 3d. The threshold value adopted should be slightly higher than the noise voltage, to avoid the influence of noise while ensuring the accuracy of instants when sensor penetrates or recedes a bubble. After that, the resulted square wave was extracted to separate the rising edge corresponding to bubble approaching a sensor and falling edge corresponding to bubble leaving a senor, as shown in Figure 3e.

The above procedures shown in Figure 3a–e are the same for *s*<sup>0</sup> to *s*3, and thus eight sets of signal are obtained. Since one bubble has two interfaces passing-by one sensor, eight rising and falling instants are produced in total (four rising edges and four falling edges), as shown in Figure 4, where the eight instants produced by one bubble are noted by *ti* and *ti* -- (*i* = 0, 1, 2, 3) with noting rising instant and - noting falling instant. To make sure each of the eight-instant group belong to the same bubble (effective bubble), the collected signals are screened by cross-checking every bubble using the method described by Equation (1):

$$\max\left( \left| t\_i'' - t\_j'' \right|\_{\mathbf{r}} \left| t\_i' - t\_j' \right| \right) < \min\left( (t\_i'' - t\_i')\_{\mathbf{r}} \left( t\_j'' - t\_j' \right) \right) \\ i, j = 0, 1, 2, 3 \text{ and i } \neq j \tag{1}$$

which guarantees the time delay when bubble approaches (or leaves) one sensor and another must be smaller than the retaining duration of one bubble. This is in the first place required by the assumption that the four-sensor probe is much smaller than the measured bubbles.

The above procedures allow omitting two kinds of ineffective bubbles. The first kind is those bubbles for which eight-instant signal is incomplete, which means that only part of the four sensors penetrate them. This means bubbles that slide or bounce off the probe. Another kind of ineffective bubbles are those whose eight-edge group does not fit the conditions of Equation (1). These bubbles can be either very small or highly deformed and are also considered as ineffective bubble. Neglecting ineffective bubbles when calculating the local flow parameters might cause error because these bubbles still have contribution to IAC. To counteract the error, the ineffective bubbles were kept counted in the program and the IAC was corrected using the average contribution of effective bubbles. The details are discussed in Section 2.4.

**Figure 4.** A group of signals produced by one bubble.

### *2.4. Deduction of Local Flow Parameters from Electrical Signals*

Before calculation of the local flow parameters, the following assumptions should be made: (1) The probe containing the four sensors is very small in size in comparison to the bubble diameter, indicating that all measured bubbles are so-called large bubbles. Small bubbles can be detectable by one or a few sensors of the four-sensor probe, but they are neglected during the signal pre-processing based on the fact that a few rising or falling edges are missing, or no edge is missing but the eight-instant does not fit Equation (1). (2) The magnitudes of interface velocity and its direction remain unchanged when an interface passes by the sensors of probe. This is true as long as the probe size is small compared to the bubbles.

The void fraction equals the ratio of the duration of the sensor contact with air to the total measurement duration. It should be noted that as there are four sensors in the probe, thus the final void fraction can be determined by their average value, as shown by Equation (2).

$$\overline{V\_f} = \frac{1}{4} \left( \sum\_{i=0,1,2,3} \frac{\sum\_{b \text{subbles}} \left( t\_i^{\prime\prime} - t\_i^{\prime} \right)}{t\_{\text{total}}} \right) \tag{2}$$

where *ttotal* denotes the total measurement time.

With the eight-instant (edge) of each effective bubble and known size and positions of the four sensors in the probe, the local flow parameters contributed by each effective bubble can be obtained. The measurement principles of other local flow parameters are schematically presented in Figure 5, which shows the relationships between different vectors.

**Figure 5.** Measurement principle of four-sensor probe. (**a**) Position of the four sensors; (**b**) Probe approaching a bubble; (**c**) Formation of measurement velocities.

Local time-averaged IAC was predicted by Ishii [57] to be related to the interface velocity projected in the normal direction of particular plane:

$$\overline{a} = \frac{1}{t\_{total}} \left( \sum\_{l} \frac{1}{|V\_l \cdot \mathbf{n}\_l|} \right) \tag{3}$$

where *<sup>l</sup>*, *V<sup>l</sup>* and *n<sup>l</sup>* denote the *<sup>l</sup>*th interface, the vector of interface velocity and the unitary vector normal to the *l*th interface, respectively, at a particular measurement point. The short line above *a* indicates a time-averaged value.

It is worth noting that, as mentioned in Section 2.3, since the neglected ineffective bubbles also contribute the IAC, Equation (3) must be corrected by multiplication factor, as shown in Equation (4):

$$
\overline{a}' = \overline{a}f = \overline{a}(\frac{N\_{\rm eff} + N\_{\rm infff}}{N\_{\rm eff}}) \tag{4}
$$

where *N*eff and *N*ineff are the number of effective and ineffective bubbles, respectively. This correction to IAC has taken the contribution of each ineffective bubbles to be equal to the average contribution of effective bubbles. A bubble is treated as ineffective mainly due to its small size and continually deforming feature. However, because of its higher surface to volume ratio, its contribution to IAC is usually larger than a bubble of large size and regular shape. Therefore, the result from Equation (4) is still expected to be lower than the true IAC.

As the methods proposed by previous researchers who have employed trigonometric functions to determine the direction of interface are not conducive to comprehensible and fast calculation, a distinct and brief vector-based calculation is proposed and performed to obtain the interface direction here in this section. The cross-product of two vectors is also a vector and its direction is perpendicular to the plane formed by the original two vectors, hence the normal vector of one interface can be determined by three velocities measurable by the probe, as shown in Figure 5c and mathematically by Equation (5):

$$m\_l = \frac{(V\_2 - V\_1) \times (V\_2 - V\_3)}{\left| (V\_2 - V\_1) \times (V\_2 - V\_3) \right|} \tag{5}$$

where *V*1, *V*<sup>2</sup> and *V*<sup>3</sup> are the three measured velocities of a particular interface, respectively. They, respectively, have the same directions with the three position vectors formed by the probe, *S*1, *S*<sup>2</sup> and

*S*3, and their magnitudes are obtained by dividing each position vector by the time delay from the probe signals. Mathematically, it can be expressed by:

$$\mathcal{V}\_{i} = \frac{\mathbb{S}\_{i}}{(t\_{i} - t\_{0})}, \quad i = 1, \ 2, \ 3 \tag{6}$$

where *ti* (*i* = 1, 2, 3) and *t*<sup>0</sup> denote the time-instant of rising or falling edges of *si* and *s*0, respectively. (*ti* <sup>−</sup> *<sup>t</sup>*0) are the time-instant delays between *si* and *<sup>s</sup>*0. When *ti* and *<sup>t</sup>*<sup>0</sup> correspond to rising edges, *<sup>n</sup><sup>l</sup>* and *V<sup>i</sup>* are the unitary normal vector and measured velocities of the front interface of a bubble, respectively; and, when *ti* and *<sup>t</sup>*<sup>0</sup> corresponds to falling edges, *n<sup>l</sup>* and *V<sup>i</sup>* are the unitary normal vector and measured velocities of the rear interface of a bubble, respectively.

The interfacial measurement theorem proposed by Shen [36] indicates the projections on the **n***<sup>l</sup>* of the interface displacement velocity vector *V<sup>l</sup>* and the measured velocities *V<sup>i</sup>* are the same, and the theorem can be expressed as:

$$\mathbf{V}\_{l'}\mathbf{u}\_{l} = \mathbf{V}\_{l'}\mathbf{u}\_{l}, \ i = 1,2,3 \tag{7}$$

from which the velocity component into the normal direction of interface can be easily obtained. In contrast, obtaining the three components of *V<sup>l</sup>* requires more assumptions or measuring parameters, for instance a further assumption of sphere-shape bubble [51] or symmetric bubble [56]. However, the application of these assumptions is only suitable when bubbles encountered in two-phase are not highly distorted or deforming.

A new method to get the whole components of *V<sup>l</sup>* is proposed in this study by using the known or legitimately assumed velocity direction, i.e., the unit vector of *Vl*. Although the shape and size of bubbles keep changing in two-phase flow, the direction of bubbles velocity is usually constant and thus the direction of *V<sup>l</sup>* often remains constant or constant on an averaged level at one fixed position of the flow field. As a result, by a known or assumed interface velocity direction *n<sup>v</sup>* (Figure 5), the magnitude of *Vl* can be obtained as follows:

$$V\_{l'} \mathfrak{n}\_{\upsilon} \cdot \mathfrak{n}\_l = V\_{l'} \mathfrak{n}\_l, \ \ i = 1, 2, 3 \tag{8}$$

The above equation applies well when the two-phase flow is limited internal flow, where the direction of the interface velocity can be regarded as parallel to the channel axis. For arbitrary multi-dimensional two-phase flow, the equation and the resulted method still apply, as long as the local flow direction can be provided by legitimate assumption, flow simulation or prior measurement. For a fixed location, the local time-averaged magnitude of interface velocity can be obtained by averaging through numerous bubbles and denoted by *Vl*.

As shown by Equation (9), the local time-averaged chord length of a bubble at a fixed position is an averaged product of the *Vl* and the averaged time duration of the four signals at high voltage.

$$\overline{C} = \frac{1}{N} \sum\_{N} \overline{V\_{I}} \frac{\sum\_{i=0}^{3} \left(t\_{i}^{\prime\prime} - t\_{i}^{\prime}\right)}{4} \tag{9}$$

where *N* denotes the number of bubbles during the measurement duration.

The direction of interface **n***<sup>l</sup>* can be shown more explicitly by the angle between **n***<sup>l</sup>* and the axis of the flow by Equation (10).

$$\theta = \cos^{-1} \frac{\mathfrak{n}\_l \cdot \mathfrak{n}\_{\text{axis}}}{|\mathfrak{n}\_l \cdot \mathfrak{n}\_{\text{axis}}|} \tag{10}$$

where *naxis* is the direction of the channel axis.

### *2.5. Innovations of the Present Probe Algorithm*

In the pioneering literature related to four-sensor probe, the interfacial velocity can only be obtained for the component that is vertical to the interface itself. If the full interfacial velocity vector is intended to be obtained, assumptions must be made. For instance, bubbles formed in the flow field in [51] were so small that they were assumed to be spherical. In [53], the authors made an assumption that the interface is very large and it moves only in its normal direction. In [56], the flattened bubbles were regarded as symmetrical to a center plane. Although these assumptions help ascertain the full components of the interfacial velocity, they only apply in special cases since the interface in practical gas–liquid flow is usually quite complicated.

Considering the fact that the explicit expression of bubble or interfacial velocity for four-sensor probe has not been completely developed, this paper proposes that, if the moving direction of interface can be ascertained prior to probe measurement, then the above assumptions will be unnecessary and the full interfacial velocity can be acquired. Fortunately, the fields of averaged moving directions of the interface, for so many quasi-steady two-phase flows, are actually the flow fields for the two-phase mixture and are easy to make certain through methods of prior measurement, simulation or even legitimate assumptions. These are the primary innovations of the present algorithm for four-sensor probe, as shown in Section 2.4.

Besides, the signal deduction process is all vector-based, which is different from the complex matrix, trigonometric functions and tensors. Although the basic rules are the same in essence and the results are expected to be same, no matter vector-based or trigonometric function-based algorithms are used, the method proposed in this paper has the merits of easy to read, clear and efficient to modify or improve.

### **3. Experimental Facilities**

To validate the availability and correctness of the probe measurement methods, an experiment was performed to compare the local flow parameters obtained by the four-sensor probe and visual measurement, including void fraction, IAC, velocity, chord length and interface direction. An air–water two-phase flow with bubbles approximately the same size injected vertically up in a transparent tube was adopted. There is no doubt that such flow properties can be easily determined by senor probe or visual techniques with high accuracy. Therefore, it was chosen as a validation case for the proposed four-sensor probe measurement method.

Although the proposed probe and algorithm have not yet been validated, the existing fundamental measurement principles described in Section 2 can also apply to micro and conventional large systems. For fierce two-phase flow scenarios, corrosive fluid, high pressure system, high flow rate and cases of flow experiencing heat and mass transfer, the measurement methods still apply as long as the probe is fixed strongly and prevented from damage. The experiment in this study was only designed for the purpose of primary method validation, considering that it is convenient to be measured by visual technique and easy to be replicated.

### *3.1. Bubbly Flow in Vertical Pipe*

An air–water two-phase bubbly flow system in a vertical tube was built to get the local flow parameter by both four-sensor probe and visual techniques. The test facility and flow rate were chosen for obtaining a simple and steady two-phase flow pattern. As shown in Figure 6, a transparent glass tube with the length of 0.5 m and internal diameter of 8 mm was used as the test section, and water was filled up to a height of 0.4 m during the experiment. Air produced by air compressor and regulated by surge tank and control valve was injected from the bottom of the tube. The air flow rate was maintained at 0.1 L/min during the test. As a result, it was found that a steady series of cap bubbles was produced inside the tube.

The top of the transparent tube was open to the atmosphere, and the four-sensor probe was vertically mounted in a one-dimensional sliding module with its tip pointing downside, so the probe could move horizontally to measure parameters across the tube diameter transversely. Only the 4.5 mm in the middle part of the 8-mm-diameter tube was accessible for measurement, resulting in a range from −2.5 to 2 mm with the interval of 0.5 mm.

The DC power supply of the probe was from a 9-V battery to avoid voltage fluctuation characterizing AC power supplier. The probe was connected to a data acquisition card which then transferred the collected data to a laptop. The data collection system was able to collect and transform the analog electrical signal to digital form at frequency of 10 kHz for each of the four sensors. For every transverse position of the probe in the tube, the data collecting persisted for 80 s and thus 800,000 data points were obtained for each sensor.

A high-speed camera was employed to record the images of bubble in a system without a probe. Images of 1262-pixel vertically and 710-pixel horizontally with a frequency of 50 fps were shot for 80 s. Through the image processing and analysis, the local flow parameters of the two-phase flow could be obtained, which is typical for the so-called visual measurement.

**Figure 6.** Experiment loop of the present study. (**a**) Schematic diagram of the experimental loop; (**b**) Real photo of the experiment loop.

### *3.2. Visual Measurement Techniques*

Besides the probe measurement described in Section 2, the flow parameters can also be obtained by the visual measurement which contains image recording and processing.

A series of continuous captured images is shown in Figure 7, with time intervals between each of 20 ms. It can be seen that the recorded bubbles are roughly in cap shape, and the space intervals between bubbles are roughly constant. When a bubble reached the location where the four-sensor probe was located, as shown by the red frame of Figure 7, the image was taken as one of the images constituting the visual measurement.

**Figure 7.** A series of images showing continuous captured bubbles.

For each of the chosen images, the subsequent bubble image processing is shown in Figures 8 and 9. For the first step of image processing (Figure 8a as an example), each bubble was cropped out, according to marked edges and converted into binary black and white image (matrix), as shown in Figure 8b–d, respectively. Five successive chosen bubbles are shown in Figure 9a, and it can be seen that, although the flow conditions remain unchanged during the test, the shape of cap bubbles change continuously. To obtain the interface direction and bubble chord length, the images of bubble were added up in MATLAB code and divided by the number of bubbles to get the averaged bubble shape. The resulted image (matrix) shown in Figure 9b stands for the probability of a pixel occupied by gas phase (void fraction). As the number of images *n* increases, the difference of the resulted images reduces to minor, and it was found that *n* = 100 is enough in this research. By binarization of the last image of Figure 9b taken with 0.5 as the threshold value, the averaged bubble and its edge are shown in Figure 9c,d, respectively. Based on Figure 9d, the time-averaged bubble interface direction θ and bubble chord length *C* at different radial location can be read by MATLAB.

**Figure 8.** Cut out bubbles and its binarization. (**a**) Bubble series; (**b**) Bubble isolation; (**c**) Bubble cropped out; (**d**) Bubble binarization.

**Figure 9.** Averaging of bubbles. (**a**) Five bubbles in a sequence; (**b**) Adding different number of bubbles; (**c**) Result of the averaged bubble; (**d**) Outer shape of the averaged bubble.

Assuming the averaged bubble in Figure 9c is axisymmetric, the bubble volume can be obtained by integrating the bubble's cross area at each horizontal layer throughout the bottom to top of the bubble, as shown in Equation (11). The total area of the interface can be obtained by integrating the bubble's interface area at each horizontal layer throughout the bottom to top of the bubble, as shown in Equation (12):

$$\overline{B\_V} = \int\_{bottom}^{top} \pi r(y)^2 dy \tag{11}$$

$$\overline{B\_{\rm I}} = \int\_{bottom}^{top} 2\pi r(y) dy \tag{12}$$

where *r*, a function of vertical location, is the radius of bubble at each horizontal layer.

The time-averaged void fraction across the whole diameter can be obtained through dividing *BV* by the average interval volume between two successive bubbles, as shown in Equation (13):

$$
\overline{\dot{V}}\_f = \frac{\overline{B\_V}}{A\overline{L}}\tag{13}
$$

where the two short lines above *Vf* indicate time-averaged and space-averaged value for the same time, A is the cross area of the test tube and *L* is the average distance between two successive bubbles.

In a similar manner, the time-averaged IAC across the whole diameter can be obtained through dividing *Ba* by the average interval volume between two successive bubbles, as shown in Equation (14):

$$
\overline{\overline{a}} = \frac{\overline{B\_a}}{A\overline{L}}\tag{14}
$$

The time-averaged bubble velocity equals the bubble production frequency multiplied by the averaged bubble distance and can be expressed by:

$$
\overline{V\_l} = \frac{N}{t\_{total}} \overline{L} \tag{15}
$$

where *N* and *ttotal* are the total number of bubbles produced and the total measurement time, respectively.

#### **4. Results and Discussion**

It is worth noting that the four-sensor probe can only give out local flow parameters of two-phase flow. It is unable to discriminate different two-phase flow regimes of stratified, slug, wavy, etc. Extra correlation research is required to make connections between local flow parameters to global flow regimes. Measuring local flow parameters at multiple locations is the purpose of the four-sensor probe and the validations against visual measurement are hence extended below.

For cap bubble occupying almost the tube diameter from −4 to 4 mm, ten radial locations ranging from −2.5 to 2 mm are measured by the four-sensor probe. The numbers of total detectable bubbles and the numbers and ratios of the effective bubbles in the 80-s measurement duration at each radial location are given in Table 1. It can be seen that 350 bubbles on average are produced in 80 s, resulting in a bubble frequency of 4.375 per second. The uncertainty of the counted bubble number changes between −5 to 6, indicating that the steady and uniform features of the bubbly flow, which is necessary for this verification test.

**Table 1.** Number and ratio of the effective bubbles to total detected bubbles.


As expected, the effective number of bubbles recognized by the MATLAB code is lower than their total number. Meanwhile, both the effective number and the ratio of the effective number to the total number decrease, and this decrease becomes faster toward the ends of the test range. For the test location beyond 2 mm, the ratio can be well below 0.5. This is because, as the probe moves away from tube axis, the bubble interfaces become more inclined and it becomes easy for the bubble to slip away from the probe. It should be noted that the errors of the obtained local flow parameters increase with the decrease of ratio of number of effective bubbles to their total number.

The void fraction from four-sensor probe by Equation (2) and visual measurement by Equation (13) are compared in Figure 10. Since the void fraction from Equation (13) is an average void fraction across the whole tube diameter, it appears as a horizontal straight line in Figure 10. In view of the bubble shape shown in Figure 9c, the bubbly flow should produce a void fraction distribution with a peak value in the tube axis and decrease towards the tube wall. This trend was successfully reflected by the probe measurement, as shown by the black square dots in Figure 10. Its arithmetic mean value shown by the dashed line agrees well with the visual measurement, with an overestimation of 8.4%.

**Figure 10.** Void fraction from probe and visual measurement.

The IAC from four-sensor probe by Equation (4) and visual measurement by Equation (14) are compared in Figure 11. Since the IAC from Equation (14) is an average void fraction across the whole tube diameter, it appears as a horizontal straight line in Figure 11. In view of the bubble shape shown in Figure 9c, the normal vector of bubble interface changes from vertical near the tube axis to pointing right-upwards near the tube wall, thus the angle between bubble velocity which is parallel to the tube and the normal vector of interface increase and the denominator of Equation (3) decreases. As a result, a larger IAC should be observed near the tube wall, with a minimum value emerging at the tube axis. This trend was successfully captured by the probe measurement, as shown by the black square dots in Figure 11. Its arithmetic mean value shown by the dashed line agrees well with the visual measurement, with an underestimation of 1.7%.

**Figure 11.** IAC from probe and visual measurement.

The bubble velocity from four-sensor probe by Equation (8) and visual measurement by Equation (15) are compared in Figure 12. Because bubble moves as an entity, it is worth noting that the bubble velocity measured at different radial locations should remain constant when the probe moves from tube axis to tube wall. However, the black square dots in Figure 12 representing probe measured velocity demonstrate fluctuation feature. This is mainly caused by error and it can be seen

the error increases towards the tube wall. Nevertheless, its arithmetic mean value shown by the dashed line agrees well with the visual measurement, with an overestimation of 9.3%.

**Figure 12.** Bubble velocity from probe and visual measurement.

The chord length from four-sensor probe by Equation (9) and visual measurement are compared in Figure 13. Both the probe measurement and visual measurement show the same trend with the chord length emerging in the middle with maximum value and decreasing towards both sides. The largest deviation between them is 8.7% at −2.5 mm, and the deviation decreases towards the pipe axis.

**Figure 13.** Bubble chord length from probe and visual measurement.

Although the interface direction is required to calculate the interface displacement velocity in the interfacial measurement theorem expressed by Equation (7), it is rarely verified against other measurement techniques in the accessible literature. Instead of directly using of the interface normal vector, the interface direction can be represented more conveniently by the angle between the interface normal vector and the tube axis. The angle from four-sensor probe by Equation (10) and visual measurement are compared, as shown in Figure 14. Figure 14a,b shows the comparisons at front interface and rear interface, respectively. As can be seen, both measurement methods show very close results with the deviation between them growing distinct towards the tube wall, and the largest errors are 22.4% and 3.1% for the front and rear interfaces, respectively.

The quantitative comparisons of the measured parameters, the deviation of the probe measurement from the visual measurement are shown in Tables 2 and 3. From the above, it is demonstrated that the measurement of bubbly flow by the four-sensor probe can give agreeable flow parameters with visual measurement techniques in aspects of void fraction, IAC, bubble velocity, bubble chord length and interface direction. Thereby, it also proves the correctness of the proposed method.

**Figure 14.** Angle between interface normal and tube axis by probe and visual measurement. (**a**) Bottom surface direction of the bubble; (**b**) Top surface direction of the bubble.



It is common knowledge that the bubble condensation is a kind of typical enhanced heat and mass transfer method encountered in numerous industrial processes such as steam–air mixture injected into subcooled zone, subcooled flow boiling and direct contact condensation [27]. Further knowledge of the bubble interface phenomena and the accurate measurement of the local flow parameters at the multi-scale interfaces, such as the void fraction, IAC, bubble velocity, bubble chord length and interface direction, are important parameters to study the process and mechanism of heat and mass transfer in gas–liquid two-phase flow. As previously reported in the literature [27,45,46,58], during the gas–liquid two-phase heat and mass transfer, the interface structure, void fraction, IAC, bubble equivalent diameter, their velocity, etc. show nonlinear variations, resulting in that they are difficult to be measured quickly and accurately. Although the double-sensor probe with multiple assumptions can also be employed to measure these parameters, it is hard to obtain realistic local data in multi-dimensional two-phase flow. As a result, the method proposed in the present study is expected to promote the solution of this problem.

**Table 3.** Probe and visual measurement of bubble chord length and interface direction.


The interface velocity and its direction obtained using the brief vector-based calculation method proposed in this study are beneficial to further master the heat and mass transfer process at the interface from a macro viewpoint. Meanwhile, the IAC is an important parameter to characterizing the interface transfer phenomenon (heat and mass transfer). Compared with the double-sensor probe, the four-sensor probe is used to measure the IAC to be able to solve the detection error caused by the escape and retreat of bubbles. On this basis, by combining the proposed vector-based signal processing method, a variety of assumptions can be avoided, which will be of great significance to further study the internal relationship between IAC and bubble size, bubble deformation and condensation effect, as well as to modify the calculation model of IAC. This shows that the proposed method can be widely applied in the future research of heat and mass transfer characteristics.

Combined with several groups of four-sensor probes and the new vector-based signal processing method, the variation laws of the void fraction, the bubble size and number of bubbles in the condensation area can be measured accurately, quickly and conveniently, and thereby the phase distribution characteristics in the condensation area can be known. Especially during the air–steam mixture bubble condensation, the phase distribution characteristics, the interface velocity and direction can reveal the mechanism affecting the interface heat and mass transfer. Furthermore, some small bubbles with high non-condensable gas concentration appear at the later stage of the air–steam mixture bubbles condensation [58]. According to the previous study [59], when the content of non-condensable air is constant, the volume change rate of bubble increases with the decrease of bubble diameter. This means that the process of small bubble condensation is of great significance for studying the enhancement of condensation heat transfer of air–steam mixture bubbles. Therefore, the proposed method in this study will also be promising for investigating the influence mechanism of non-condensable air at the interface on the tiny steam bubble condensation enhancement.

### **5. Conclusions**

The void fraction, IAC, bubble size, velocity, etc. of two-phase flow have great influences on the heat and mass transfer characteristics, the reaction efficiency and the operation safety of multiple chemical applications. To obtain detailed knowledge of two-phase flow, a miniaturized four-sensor probe was made firstly in this study. Then, a new method based on vector calculation to get the interface direction and a new method to get the bubble velocity magnitude were proposed, by which other local flow parameters can be obtained from the signals produced by four-sensor probe. For the purpose of verification of the probe made and the elaboration algorithm developed, an experimental facility was also built and air–water two-phase bubbly flow was tested. Besides the measurement by probe, the flow parameters were also obtained by video recording and image process techniques. For the first time, the direction of bubble interface was compared with visual image of bubble.

The calculation of interface velocity vector was shown to be realizable using legitimate interface velocity direction. The comparison between the results obtained from probe measurement and visual measurement indicates their good agreement. However, since the number of ineffective bubbles increases moving towards tube wall, the measurement error by probe also increases while the probe moves away from tube axis. The averaged values across the tube are over predicted by 8.4%, under predicted by 1.7% and over predicted by 9.3% with probe measurement for void fraction, IAC and bubble velocity, respectively, in comparison to the visual measurement counterpart. The chord length and angle between tube axis and normal to the interface show great consistency with visual results, with maximum deviations of 8.7% and 22.4%, respectively.

After more validation tests against multiple kinds of two-phase flow parameters, it can be concluded that the methods proposed in this study are promising in the characteristics and mechanism studies of the gas–liquid heat and mass transfer, the direct measurement of gas–liquid two-phase flow and providing significant database for the improvement of two-phase flow models.

**Author Contributions:** Conceptualization, X.Q. (Xiaohang Qu) and Y.Z.; methodology, X.Q. (Xiaohang Qu), Q.G. and Y.Z.; software, X.Q. (Xiaohang Qu) and X.Q. (Xiaoni Qi); validation, X.Q. (Xiaohang Qu), Q.G. and Y.Z.; formal analysis, X.Q. (Xiaoni Qi), Q.G., L.L. and Y.Z.; investigation, X.Q. (Xiaohang Qu) and Y.Z.; data curation, X.Q. (Xiaohang Qu), Y.Z. and L.L.; writing—original draft preparation, X.Q. (Xiaohang Qu), and Y.Z.; writing—review and editing, Q.G. and X.Q. (Xiaoni Qi); visualization, X.Q. (Xiaohang Qu), Y.Z., X.Q. (Xiaoni Qi) and L.L.; and supervision, Q.G. and Y.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant Nos. 51806128 and 51879154) and the Natural Science Foundation of Shandong Province (Grant Nos. ZR2019BEE008 and ZR2018LE011).

**Acknowledgments:** The authors gratefully acknowledge financial support from the National Natural Science Foundation of China (Grant Nos. 51806128 and 51879154) and the Natural Science Foundation of Shandong Province (Grant Nos. ZR2019BEE008 and ZR2018LE011).

**Conflicts of Interest:** The authors declared that there is no conflict of interest.

### **Nomenclature**

The following abbreviations are used in this manuscript:


### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **CFD Analysis of Subcooled Flow Boiling in 4** × **4 Rod Bundle**

### **Ye-Bon Seo, SalaiSargunan S Paramanantham, Jin-Yeong Bak, Byongjo Yun and Warn-Gyu Park \***

School of Mechanical Engineering, Pusan National University, Busan 46241, Korea; seoyebon308@pusan.ac.kr (Y.-B.S.); salaisargunansp@pusan.ac.kr (S.S.P.); jinyb0309@pusan.ac.kr (J.-Y.B.); bjyun@pusan.ac.kr (B.Y.)

**\*** Correspondence: wgpark@pusan.ac.kr

Received: 10 June 2020; Accepted: 24 June 2020; Published: 30 June 2020

**Abstract:** Rod bundle flow is an important research field related to reactor cooling in nuclear power plants. Owing to the rapid development of computerized performance assessments, interest in coolant flow analysis using computational fluid dynamics has garnered research interest. Rod bundle flow research data compared with experimental results under various conditions are thus needed. To address this, a boiling model verification study was conducted with reference to experiments. This study adopts the Reynolds-averaged Navier–Stokes equation, a practical analysis method compared to direct numerical simulation and large eddy simulation, including turbulence modeling, to predict the flow of coolant inside a rod bundle. This study also investigates void behavior in low-pressure subcooled flow boiling using a Eulerian approach (two-fluid model). The rod bundle has a length of 0.59 m and a hydraulic diameter of approximately 14.01 mm. At the cross-section at a height of 0.58 m, near the exit, numerical results were compared with the experimental values of the volume fraction of vapor and interfacial area concentration. The simulation results showed good agreement with the experimental data for six different initial conditions with constant density.

**Keywords:** subcooled flow boiling; two-fluid model; vapor volume fraction; rod bundle; STAR-CCM+

### **1. Introduction**

Computational fluid dynamics (CFD) is being increasingly used in the analysis of three-dimensional two-phase flow for the safe operation of a nuclear power plant. In recent years, subcooled boiling in the rod bundle is one of the important phenomena observed in small break loss of coolant accidents (SBLOCA) for pressurized water reactor (PWR) and is also expected in the loss of heat sink accident of nuclear spent fuel pool (SFP) as a safety issue. Thus, experiments and measurements on multiphase turbulent flow in pipes or rod bundles have been conducted, and numerous fluid properties have been studied. In recent years, various simulations have been conducted in various papers [1–5]. Numerical analysis of subcooled two-phase flow at high pressure has confirmed that the Eulerian multiphase model including a Reynolds stress turbulence model is suitable for predicting the actual phenomenon [6–10]. Single-phase convective heat transfer and vapor–water two-phase boiling heat transfer experiments have been conducted in a vertical seven-rod bundle at a low flow rate near atmospheric pressure. In this case, bubble generation is dominantly determined by the enhancement factor of forced convection by the inlet mass flux and the suppression factor of nucleate boiling. Thus, it was necessary to react sensitively to the fluid velocity, and in the next step, excessive bubble prediction is prevented through a nuclear boiling suppression factor. By taking in account the enhancement factor of forced convection and the suppression factor of nucleate boiling, a new correlation in terms of boiling number, Reynolds number, and Martinelli number was found to be effective in predicting two-phase-flow boiling heat transfer coefficients [11]. Analytical studies using various turbulence models for rod bundles in various shapes have also been conducted. A 3 × 3 rod bundle flow was

analyzed using the k–ε model, shear stress transport k-ω model, and Reynolds stress model; and the pressure rate, flow velocity, and temperature were predicted [12]. For supercritical pressure conditions, the importance of the turbulence mixing factor was studied thermal hydraulics characteristics of 2 × 2 rod bundle by numerical analyzing [13]. In total of 2000 mm length of the 4 × 4 rod bundle with a mixing grid installed at 2/3 height, turbulence was processed in a Reynolds-averaged Navier–Stokes method using a nonlinear and viscous model to simulate CFD [14]. A study on the effect of spacer grid distance considered the influence of spacer grid blockage ratio and spacer grid spacing on heat transfer in a vertical rod bundle [15]. Reynolds number of the working fluid in this paper is 24,800 to 30,100 to investigate the flow effect on heat transfer. As such, this study was carried out with reference to other studies on the two-phase model, turbulence model, and various rod bundle shapes and grids.

The present study addresses the local flow inside the rod bundle. As the inside of a bundle of nuclear fuel rods has a very narrow flow path, a two-phase flow study using three-dimensional heat transfer analysis and a boiling model is required for reactor safety design and analysis. For two-phase flow analysis, the bubble size is a very important factor for accurately predicting bubble behavior in the flow path. In line with this focus, this study analyzes the results of the numerical analysis of the volume fraction of vapor (VF) and interfacial area concentration (IAC) compared with experimental values. VF is the volume of vapor in volume rate and IAC is defined by VF and interaction length scale. It focuses on the analysis of coolant boiling in a narrow flow path inside a nuclear fuel reactor 4 × 4 rod bundle. This phenomenon is solved by the Eulerian multiphase model as a segregated flow. The important elements to accurately predict the behavior of vapor are the interfacial momentum transfer force and wall boiling model of heat and mass transfer force. The usability of the detailed model used in this study has been evaluated by researchers as follows.

Tomiyama (1998) studied simple but reliable correlations for a drag coefficient by comparing experimental data. The correlations were obtained for the drag coefficient of single bubbles under a wide range of fluid properties, bubble diameters, and acceleration of gravity. These were developed according to the balance of forces acting on the bubbles in stagnant liquid, and the available empirical correlations of the terminal rising velocities of single bubbles [16].

It has been shown that in many flows, the inertial force can simply be added to the generalization of the lift force introduced by Auton (1987) [17]. The lift force of Auton (1988) consists of the rotational force and inertial force or additional mass [18].

Owing to the neglect of the wall area force, the predicted pore profile peaked instead of approaching zero at the wall. Therefore, Achard & Cartellier (1985) concluded that a different kind of transverse force is needed to break the bubbles away from the wall [19]. Later, Antal (1991) developed a two-fluid model of multidimensional laminar bubble two-phase flow to analyze vertical pipe flow. In this study, the wall lift force for the dispersed phase (i.e., bubbles) and the wall lubrication force about the wall force were calculated. The results showed good agreement with the data when the considered model was used [20].

Lahey et al. (1993) used *CTD* = 1 to predict the lateral phase distribution of a three-dimensional water/air bubble upflow or downflow in a pipe using the k–ε turbulence model and an algebraic stress law. It can be confirmed [21]. These results are very interesting because a properly formulated multidimensional two-fluid model has the inherent ability to predict the lateral phase distribution phenomena in simple and complex geometric conduits for bubbles of size 1 mm < *Db* < 6 mm.

The model of Kurul and Podowski (1990) was used to determine the influence of bubble wall area fraction properties [22]. This model can estimate the fraction of the wall area affected by the sweeping of the liquid inflow under the bubble departure.

In the two-fluid model, active nucleation site density is an important factor for estimating the IAC. Thus, the model by Hibiki–Ishii (2003) proposed that the active nucleation site density is a function of the critical cavity size and contact angle [23]. This model clearly demonstrated the dependence of active nucleation site density on the superheated wall by various studies.

Kocamustafaogullari (1983) reported that the bubble departure diameter was associated with water over a wide range of pressures [24]. Using the correlation proposed by Kocamustafaogullari, the mean deviation of the experimental data is approximately ±33%. It accurately simulated the available experimental water data, considering various free surfaces.

Cole (1960) demonstrated through experimental measurements that the main forces acting on bubbles leaving the surface at high heat flows are the buoyancy and drag forces [25]. The dimensionless relationship was developed by the bubble velocity, bubble diameter, and contact angle at breakoff.

There are many other studies. However, the boiling phenomenon of the subchannel is still a field that requires various studies depending on the total length, hydraulic diameter, flow rate, and heat transfer rate. This study analyzes the coolant flow in 4 × 4 rod bundles with applied symmetrical conditions. The diameter, pitch, and hydraulic diameter of the rods are 16, 21, and 14.01 mm respectively, and the total length is 590 mm. The cross-sectional area of the coolant is 0.004008 m2. The grid size has three prism layers total of 0.8 mm near the wall, and 1 mm in other areas. Since the fluid belongs to high-Re, High-y+ wall treatment (applicable to 30 ≤ y+ ≤ 100) was applied. This simulation is calculated with steady state and constant density. Initial condition covered a mass flux ranging from 250 to 310 kg/m2s, heat flux ranging from 148 to 200 kW/m2, and inlet subcooled temperature ranging from 20.2 ◦C to 25.2 ◦C. Numerical analysis is performed by STAR-CCM+ ver. 12.06.011 considering three-dimensional turbulence.

### **2. Subcooled Boiling Experiment in 4** × **4 Rod Bundle**

In this study, an experiment was conducted in the subcooled boiling water flow with a 4 × 4 rod bundle heater. Instruments such as Coriolis flow meters, static pressure gauges, differential pressure gauges, thermal couples, and power meters were equipped in the facility, and control valves for precisely controlling the flow rate and preheater to maintain the target inlet temperature were installed. The measurement instrument used in the experiment and the uncertainty are listed in Table 1. The pressure could be adjusted with a pressurizer connected to the water storage tank consisting of a test loop. A separator for separating steam and water was installed at the outlet of the test section, and a condenser and heat exchanger for heat removal from steam and water, respectively, were installed at the rear end of the water separator. An additional heat exchanger was installed at the outlet of the test section to prevent the flashing of high-temperature water. To minimize heat loss from the flow test channel during the experiment, insulating material was applied to the outside wall of the flow channel. The test section was designed as a square channel of which the length of the sides was 85 mm. Each heater rod was 16 mm in outer diameter and 590 mm in heated length (Zh). The pitch-to-diameter ratio was 1.3 in the bundle arrangement; therefore, the hydraulic diameter (DH) of the whole subchannel was 14 mm.


A two-sensor optical fiber probe (2S-OFP) consisting of one front and one rear sensor was applied to measure local two-phase flow parameters in the rod bundle. The local two-phase flow parameters measured in the experiment were the void fraction (α*<sup>t</sup>* ), interfacial area concentration (*ai t* , IAC), interfacial velocity (*vij t* ), and bubble diameter (*Dsm t* ). The local void fraction is the ratio of passage time of bubbles per unit time at a local point in the two-phase flow condition. The IAC, defined as the interfacial area per unit volume, is a parameter of prime importance in terms of mass, momentum, and energy transfers between phases. In this study, the interfacial velocity is the velocity of the actual

velocity of the bubbles along the axial direction of a flowing channel. The detailed equations for each parameter are as follows.

$$\overline{\alpha}^{\mathbf{t}} = \frac{1}{\mathbf{T}} \sum\_{\mathbf{j}}^{\text{N}\_{\text{b}}} (\mathbf{t}\_{\text{F2}} - \mathbf{t}\_{\text{F1}}) \tag{1}$$

$$\overline{\mathbf{a}\_{\mathrm{i}}}^{\mathrm{t}} = 2\mathbf{N}\_{\mathrm{b}} \xrightarrow[\overline{\mathrm{V}\_{\mathrm{ij}}}]{} \mathrm{I}\Big(\boldsymbol{\phi}\_{\mathrm{j,max}}\big)\tag{2}$$

$$\overline{\mathbf{v}\_{\text{ij}}}^{\mathbf{t}} = \frac{\mathbf{N}\_{\text{b}}}{\sum\_{\text{j}}^{\mathbf{N}\_{\text{b}}} \frac{(\mathbf{t}\_{\text{R}1} - \mathbf{t}\_{\text{F1}})}{\mathbf{S}}} \tag{3}$$

where T is the total measurement time, Nb is the number of bubbles in contact with local optical fiber sensor, tF2 is the time when the bubble completely passes through the front sensor, tF1 is the time when the bubble starts to pass through the front sensor, tR1 is the time when the bubble starts to pass through the rear sensor, I <sup>φ</sup>j,max is the angle between the mean flow direction vector and the velocity vector of the interface, proposed by Hibiki et al. (1998) [26], and superscript t is the time average. The bubble diameter is six times the void fraction over the interfacial area concentration. Shen and Nakamura (2013) [27] showed a measurement error of ±7.8% when measuring the void fraction of small bubbles with an optical fiber probe technique. The interfacial velocity for small bubbles estimated with error less than ±2.3% when comparing the photographic and double sensor probe method used in our laboratory.

The measurement was performed on a total of 26 points in the area corresponding to 1/8 of the center subchannel at Zh/DH = 41.8 downstream from the inlet of the test section. The flow condition covered a mass flux ranging from 250 to 310 kg/m2s, heat flux ranging from 148 to 200 kW/m2, and an inlet subcooled temperature ranging from 20.2 ◦C to 25.2 ◦C at a 200 kPa inlet pressure.

### **3. Numerical Methods**

For the analysis performed in this study, a two-fluid analysis method based on the Eulerian equation and a standard k–ε turbulence model were used [28,29]. In particular, this study considers the drag force of Tomiyama et al. (1998) [16], the lift force of Auton et al. (1988) [18], the wall lubrication force of Antal et al. (1991) [20], the turbulent dispersion force of Lahey et al. (1993) [21], and the virtual mass force of Auton et al. (1988). [18] The wall heat distribution model of Kurul and Podowski (1990) [22] was used to calculate the boiling flow, considering the interfacial length scale from Bak (2015) [30] of the bubble diameter model. Hibiki–Ishii (2003) [23] nucleation density, Kocamustafaogullari (1983) [24] bubble departure diameter, and the Cole (1960) [25] bubble departure frequency model were used to calculate the wall bubble generation.

### *3.1. Turbulence Model*

Prior to the numerical analysis, the standard k–ε and realized k–ε turbulence models were compared to confirm the turbulence model sensitivity [31]. Figure 1 shows that there is no significant difference between the results of the two models. However, when calculations were made using the standard k–ε model under the same conditions, the results were closer to the experimental values. Therefore, the standard k–ε model was adopted for all cases.

**Figure 1.** Turbulence model test.

### *3.2. Two-Phase Interaction Models*

The drag force and coefficient of Tomiyama et al. (1998) [16] are given by

$$F\_{ij}^D = \mathbb{C}\_{ij}^D \frac{1}{2} \rho\_c |u\_j - u\_i| (u\_j - u\_i) \left(\frac{a\_{ij}}{4}\right) \tag{4}$$

$$\mathcal{L}\_{ij}^{D} = \max \left[ \left( \frac{24}{Re} (1 + 0.15 Re^{0.687}) \right), \frac{8Eo}{3(E0 + 4)} \right] \tag{5}$$

In Equation (5), *Re* is the Reynolds number and *Eo* is the Eotvos number, expressed as

$$\mathcal{R}\mathfrak{c} = \frac{\rho\_{\mathfrak{c}}|\mu\_{\mathfrak{r}}|l}{\mu\_{\mathfrak{c}}} \tag{6}$$

$$Eo = \frac{\left|\rho\_c - \rho\_d\right| \mathcal{S} l^2}{\sigma} \tag{7}$$

This method is suitable for liquid-gas bubble systems with a low Reynolds number and low Morton number flow; for example, small bubbles in an air–water system. It is only available when the dispersed phase is a gas. The lift force and coefficient of Auton et al. (1988) [18] are given by

$$F\_{ij}^{l} = \mathbb{C}\_{ij}^{l} \rho\_{\mathbb{C}} \alpha\_{d}[\mu\_{r} \times (\nabla \times \mu\_{r})] \tag{8}$$

$$\mathcal{C}\_{ij}^{l} = \text{Const.} (-0.5 \sim 0.5) \tag{9}$$

The wall lubrication force that prevents bubbles from touching the wall, and the coefficient of Antal et al. (1991) [20] are given by

$$F\_{ij}^{\rm NL} = -\mathcal{C}\_{ij}^{\rm NL} \frac{\alpha\_d \rho\_c u\_r^2}{d} n \tag{10}$$

$$\mathbf{C}\_{ij}^{\text{NL}} = \max \left[ \mathbf{C}\_1 + \mathbf{C}\_2 \frac{d}{y\_w}, 0 \right] \tag{11}$$

In Equation (11), *C*<sup>1</sup> and *C*<sup>2</sup> are constant. The turbulence dispersion force and coefficient of Lahey et al. (1993) [21] are given by

$$F\_{ij}^{TD} = \mathbb{C}\_{ij}^{TD} \frac{1}{2} \rho\_i (u\_j - u\_i) \Big(\frac{a\_{ij}}{4}\Big) \frac{v\_c^t}{\sigma} \Big(\frac{\nabla a\_j}{a\_j} - \frac{\nabla a\_i}{a\_i}\Big) \tag{12}$$

$$\mathcal{L}\_{ij}^{TD} = 1 \tag{13}$$

The virtual mass force and coefficient of Auton et al. (1988) [18] are given by

$$F\_{ij}^{VM} = \mathbb{C}\_{ij}^{VM} \rho\_{\mathbb{C}} \alpha\_d (a\_j - a\_i) \tag{14}$$

$$C\_{ij}^{VM} = 0.5\tag{15}$$

The nucleation site number density determines the number of locations on the heated surface where bubbles form, per unit area. The nucleation site number density of Hibiki–Ishii (2003) [23] was used as follows

$$d\_{\omega} = d\_1 \theta \left(\frac{\sigma}{g \Delta \rho}\right)^{0.5} \left(\frac{\sigma}{\rho\_{\mathcal{S}}}\right)^{0.9} \tag{16}$$

where *<sup>d</sup>*<sup>1</sup> is 2.64 <sup>×</sup> <sup>10</sup>−<sup>5</sup> *<sup>m</sup>*/*deg*, and <sup>θ</sup> is 41.36◦ for the active nucleation site number density as a function of critical cavity size and contact angle. It is designed for use with the Kocamustafaogullari (1983) [24] model for bubble departure diameter as follows

$$n^{\prime\prime} = \overline{n^{\prime\prime}} \left( 1 - \exp\left[ -\frac{\theta^2}{8\mu^2} \right] \right) \left( \exp\left[ f(p^+) \frac{\lambda^{\prime}}{R\_c} \right] - 1 \right) \tag{17}$$

The bubble departure frequency model of Cole (1960) [25] was used as follows

$$f = \sqrt{\frac{4}{3} \frac{g(\rho\_l - \rho\_{\mathcal{S}})}{d\_{o} \rho\_l}}\tag{18}$$

### *3.3. Bubble Diameter Model*

In this three-dimensional two-phase flow analysis, the bubble size is a very important factor for accurate prediction of bubble behavior in the flow path. Yao and Morel (2004) and Yeoh and Tu (2005) developed a kinematic model based on the IAC and bubble number density transport equations for accurate bubble size prediction. The model has been developed and requires subsidiary equations [32,33].

The bubble size is known to be related to the hydraulic diameter of the flow path, Laplace length scale, volume fraction, and Reynolds number of the bubble through the previously developed models of Zeitoun et al. (1994) [34], Zeitoun and Shoukri (1996) [35], and Hibiki et al. (2006) [36]. Thus, bubble size model is formed by

$$d\_{sm} = \alpha^A \mathcal{N}\_{Reb}^B L\_O^C \tag{19}$$

where α is the volume fraction. The bubble Reynolds number (*NReb*) to consider the turbulent dissipation term (ε) and the Laplace length scale (*LO*) are given by

$$N\_{Rcb} = \frac{\varepsilon^{1/3} L \sigma^{4/3}}{\upsilon\_f} \tag{20}$$

$$L\_{\bullet} = \sqrt{\frac{\sigma}{g\Delta\rho}}\tag{21}$$

In Equation (20), the subscript f means liquid phase. In Equation (21), σ is the surface tension, *g* is the gravitational acceleration, and Δρ is the density difference between phases.

Based on Equation (19), Hibiki realized a bubble size model applicable to the bubble flow as

$$d\_{\rm sm} = 4.34 \alpha^{0.36} N\_{\rm Reb}^{-0.696} L\_{\rm O}^{0.665} / D\_{\rm H} \tag{22}$$

A simple form of Equation (20) is derived by adding the density ratio and excluding the hydraulic diameter included in the bubble Reynolds number [37]. Bak derived a bubble diameter model from Yun by applying DEDALE, Hibiki, and DEBORA experimental data for both low- and high-pressure conditions. The bubble size model of Bak (2015) [30] used in this study is given by

$$d\_{\rm sun} = 32.39 \alpha^{0.289} N\_{\rm Reb}^{-0.529} N\_{\rho}^{-0.06} L\_O \tag{23}$$

where the density ratio (*N*ρ) to reflect pressure is given by

$$N\_{\rho} = \frac{\rho\_f}{\rho\_{\S}} \tag{24}$$

and the subscripts *f* and *g* are the liquid phase and gas phase, respectively.

### **4. Operating Conditions**

The numerical analysis of the rod bundle has six initial conditions. It also proves that the grid test was conducted before showing the numerical analysis results.

### *4.1. Initial and Boundary Conditions*

The total length of the rod bundles was 590 mm. The diameter, pitch, and hydraulic diameter are 16, 21, and 14.02 mm, respectively. The cross-sectional area through which the coolant flows is approximately 0.00401 m2. The experimentally measured data on Lines 1–5 comparing the numerical analysis results at the cross-section of 580 mm in height are shown in Figure 2.

**Figure 2.** Measured points on Lines 1–5 for numerical analysis.

The initial conditions for the six experimental cases used in the calculations are listed in Table 2. At the outlet of all cases, the fluid has a pressure of 194.45 kPa. This shows a decrease of about 5 kPa from the inlet. Calculation was performed with constant density.


**Table 2.** Initial conditions of six cases.

In Table 2, the different conditions are highlighted in gray. Cases 2 and 3 have different inlet temperatures compared with Case 1. Cases 4 and 5 have different heat flux compared with Case 1. Case 6 has a different inlet velocity from that of Case 3.

### *4.2. Grid Test*

Given that the fluid does not have different results at a centerline of 4 × 4 and 2 × 2, the numerical analysis of the rod bundle was modeled by constructing a shape with one 2 × 2 side and by applying a symmetrical condition twice. The geometry part is the gray area in Figure 3b.

**Figure 3.** Shape of the rod bundle: (**a**) three-dimensional model viewed in the x–z plane and (**b**) cross-section of 0.58 m in height and wall conditions.

The analysis results are shown in Figure 4, with approximately one million, two million, and four million grids. A comparison of the analysis values was made between Line 1 and Line 4.

**Figure 4.** Comparison of one million, two million, and four million grids for tests on Lines 1 and 4: (**a**) volume fraction and (**b**) interfacial area concentration.

The grid test showed no significant difference between the approximately one million, two million, and four million grids. Therefore, numerical analysis was performed with one million grids, which had the shortest calculation time. The cross-sectional shape of the one million grid is shown in Figure 5.

**Figure 5.** One-million-grid rod bundle: (**a**) 3-D near the outlet and (**b**) cross-section of the outlet.

The one million grid consists of hexahedron grids. By applying three prism layers to the wall surface, the spacing is approximately 0.2667 mm near the wall, for a total of 0.8 mm. In other areas, the grid spacing was 1 mm. This grid method calculates three important measures of mesh quality at the start of a run: aspect ratio—maximum 0.99: good; skewness angle—maximum 4.79: good; volume change—maximum 0.98: good. Mean aspect ratio is 0.81, mean skewness is 4.75, mean volume change is 0.8. In terms of overall mesh quality minimum quality is 0.53, maximum 0.99, and mean 0.97.

### **5. Results**

### *5.1. Results*

Before the results of the six cases are shown, Figure 6 illustrates a stable flow falling by approximately 5.53 kPa of pressure from about 200 kPa from the inlet to the outlet. The pressure distribution is shown in Figure 6.

**Figure 6.** Pressure distribution of numerical analysis results in Case 1.

For the parameters of the results, the meanings of the volume fraction and IAC are as follows.

$$Volume\,\,fraction = \frac{Volume\,\,of\,\,vapor}{Volume} \tag{25}$$

$$IAC = 6 \frac{Volume \, fraction}{\text{Interaction length scale}} \tag{26}$$

In Equation (26), the interaction length scale is defined as the bubble diameter.

Figures 7–12 show the comparison of the results of Lines 1–5 in Cases 1–6. In each figure, (a) shows the volume fraction and (b) shows the IAC.

In Case 1, it is considered that the VF and IAC were accurately predicted. Cases 2 and 3 (subcooled temperature of 20.2 ◦C or 25.2 ◦C) have different inlet temperature from Case 1 (subcooled temperature of 22.7 ◦C). The axis shows numerical values similar to experimental values, but has smaller values to the rod surface. There was no difficulty in interpreting the same tendency, such as the slight decrease in numerical results near the rod surface, as the experimental results. The result of Case 4 shows accurate prediction. Case 4 also indicates low VF and IAC, the same as the experiments, because of the low heat flux of 148 kw/m2. In Case 5, which had a high heat flux of 200 kw/m2, the IAC was low near the rod surface, but it can be said that it measures some reasonable value in the VF. Lastly, Case 6, which had a low flow rate, was under-interpreted near the rod surface except on Line 1, but it shows that it interprets the appropriate values near the axis and across Line 1.

**Figure 7.** Case 1 numerical results of Lines 1–5 compared with experiment results: (**a**) volume fraction and (**b**) interfacial area concentration.

**Figure 8.** Case 2 numerical results of Lines 1–5 compared with experiment results: (**a**) volume fraction and (**b**) interfacial area concentration.

**Figure 9.** Case 3 numerical results of Lines 1–5 compared with experiment results: (**a**) volume fraction and (**b**) interfacial area concentration.

**Figure 10.** Case 4 numerical results of Lines 1–5 compared with experiment results: (**a**) volume fraction and (**b**) interfacial area concentration.

**Figure 11.** Case 5 numerical results of Lines 1–5 compared with experiment results: (**a**) volume fraction and (**b**) interfacial area concentration.

**Figure 12.** Case 6 numerical results of Lines 1–5 compared with experiment results: (**a**) volume fraction and (**b**) interfacial area concentration.

Figure 13 shows the numerical analysis results for the volume fraction distribution contour of all cases for the cross-section with a height of 580 mm. In Figure 13c,d, the volume fraction distributions are characterized by a low value between the rod surface and the wall. However, in Figure 13a,b,e,f, the volume fraction distributions are characterized by a high value between the rod surface and the wall.

**Figure 13.** *Cont*.

**Figure 13.** Volume fraction distribution on the cross-section with a height of 580 mm: (**a**) Case 1, (**b**) Case 2, (**c**) Case 3, (**d**) Case 4, (**e**) Case 5, and (**f**) Case 6.

This can be explained by the velocity and temperature of the liquid. Under the same conditions, Case 4 had the smallest heat flux (148 kW/m2) and Case 5 had the largest heat flux (200 kW/m2). The detailed analysis of this is shown in Figures 14 and 15, taking Cases 4 and 5 as examples. As shown in Figure 14a, the liquid velocity decreased more as it approached the outlet in the wall edge. Therefore, the liquid temperature increased because of heat transfer. However, because Case 4 was a case with very low heat flux, the temperature at the wall edge excluding the rod surface did not rise sufficiently to the boiling point, as shown in Figure 14b. Therefore, as shown in Figure 14c, air bubbles were not generated at the corners. Because no bubbles were generated, the velocity of the bubbles remained low, as shown in Figure 14d.

**Figure 14.** Distribution contour in Case 4: (**a**) velocity of liquid, (**b**) temperature of liquid, (**c**) volume fraction, and (**d**) velocity of vapor on the vertical plane at the wall edge.

**Figure 15.** Distribution contour in Case 5: (**a**) velocity of liquid, (**b**) temperature of liquid, (**c**) volume fraction, and (**d**) velocity of vapor on the vertical plane at the wall edge.

Figure 15a, however, shows that the liquid flowed with a low velocity at the wall edge from the inlet to the middle height of the rod bundle. Therefore, as the liquid velocity decreased, the temperature increased enough to reach the boiling point from the middle height, as shown in Figure 15b. Consequently, more bubbles were generated, as shown in Figure 15c. Figure 15d shows that the velocity of the liquid increased more as the generated vapor velocity increased.

### *5.2. Model Validation*

Various initial conditions were used for the numerical analysis of the six cases, with inlet fluid velocity of approximately 0.260417 or 0.322917 m/s, heat flux of 148–200 kw/m2, and inlet subcooled temperature of 20.2–25.2 ◦C. The graphs on the cross-section near the outlet from the rod can confirm the validity of the model in Figures 7–12. In particular, this model predicts the reasonable value on center of rod bundles in a wide flow path and the high inlet temperature or the large heat flux; the sub-cooled temperature is 20.2 ◦C or the heat flux is 200 kw/m2. All the volume fractions and IACs in Figures 7–12 tend to be similar to the experimental values, which indicates that this study will be useful for understanding the vapor behavior of the entire steady-state flow in the simulations for predicting the boiling phenomenon. The results shown in Figures 14 and 15 are expected to aid in predicting the correlation between liquid and volume fraction, which is the most important factor considered in bubbly flow CFD research. These CFD results generally reproduced the trend of experimental values. However, this simulation has the same kind of error in all cases. Thus, we must describe the error. Numerical analysis results in higher errors as the mass flux is higher and the heat flux and temperature are lower. Moreover, it can be confirmed that the VF is smaller than the experimental value. The analysis of this cause from the model aspect shows that is due to the use of a constant value for each coefficient of the interface momentum transfer force models in this paper. In each model, the coefficient was inputted as a certain and reasonable constant with reference to the papers, but it is not a value derived from the exact same experiment and the same working fluid, so an error occurs. In the interfacial momentum transfer force model that predicts the movement of bubbles, the most dominant force is the drag force, and then the lift force. When looking at the horizontal cross section of a vertical pipe, the lift force, which determines the location of the bubble, is one of the important factors in determining how well the bubble moves after generation and departure from the heating surface [38–40].

### **6. Conclusions**

In this study, numerical analysis of two-phase flow boiling by heat transfer was conducted. The experimental values and CFD results showed good agreement in all six cases considered. The analysis was a steady state of 3D and the calculated model was performed at constant density. The fluid was calculated using the Eulerian multiphase model including a Reynolds stress turbulence model. The interfacial momentum transfer force model and wall boiling model are important models in two-phase flow. The momentum transfer force was calculated by considering the drag force, lift force, wall lubrication force, turbulent dispersion force, and virtual mass force. The nucleation site number density, bubble departure diameter, and bubble departure frequency of the wall boiling model were considered in the heat and mass transfer force. This means that the model considered the behavioral characteristics of the fluid to be as realistic as possible. Thus, it is possible to sufficiently interpret the code. In particular, it shows stable analysis results in all cases with different conditions. The main contributions of this study are as follows.


**Author Contributions:** Conceptualization, W.-G.P.; formal analysis, Y.-B.S. and S.S.P.; funding acquisition, W.-G.P.; investigation, Y.-B.S. and S.S.P.; methodology, W.-G.P.; project administration, W.-G.P. and B.Y.; Supervision, W.-G.P. and B.Y.; writing—original draft, Y.-B.S. and J.-Y.B.; writing—review & editing, Y.-B.S., W.-G.P. and B.Y.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by "Human Resources Program in Energy Technology" of the Korea Institute of Energy Technology Evaluation and Planning (KETEP), granted financial resource from the Ministry of Trade, Industry & Energy, Republic of Korea. (No. 20184010201660) and the Nuclear Research & Development Program of the NRF (National Research Foundation of Korea) grant funded by the MSIT (Ministry of Science and ICT). (NRF-2019M2D2A1A03056998).

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

### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel: +41 61 683 77 34 www.mdpi.com

ISBN 978-3-0365-5470-9