CFD Modeling of Multiphase Flows in the Gas Turbine Engines Oil Cavities

The choice of the structure of the mathematical model of thermal-hydraulic processes in the oil cavities of the GTE rotor supports has been substantiated. A three-dimensional CFD model has been built to calculate multiphase currents involving information on the flow distribution and heat exchange given in the scientific literature. We have considered the approaches and individual models used for these purposes. The resulting solutions are consistent with the results of the experiment on a model support and generally accepted ideas about the processes in a given class of devices. The distribution of oil in the chamber, the phase current lines, the temperature and velocity fields have been given, as well as the velocity vectors for various CFD models (VOF, Euler, Inhomogeneous) and solver types (steady and non-steady). Based on the analysis of the results obtained, it has been found that the Euler model involving a non-steady solver yields the smallest difference with the experimental values for a heat transfer coefficient. In all cases, when gravity is considered, there is an asymmetrical distribution of the oil film. The result is a change in the thermal resistance of the boundary layer and, consequently, in the heat transfer coefficient along the bearing chamber circumference. This largely determines the heat flow through the chamber wall. The proposed method of modeling workflow in the support's oil cavity is based on a mathematical notation of the heterogeneous monodisperse oil-air flow with an algorithm of inversion of the structure of two-phase flow in the near-wall region from the drip into the bubble. That makes it possible to more accurately calculate the temperature states of the GTE rotor support elements and the system that ensures the proper operation of the bearing by correctly determining the heat transfer coefficient on the part of the oil-air mixture. The constructed model makes it possible to numerically investigate the applicability of those known and to derive the new correlation dependences for the mean value of aheat transfer coefficientin the oil cavity of rotor support that is used in engineering calculations. The model also makes it possible to numerically investigate the impact of the geometry, the rotor rotation frequency, and the phases flow rates on heat output in the oil cavity.

most complex and little-studied ones. A proper description of processes in the support's oil cavity, also referred to as a "bearing chamber," is the determining factor in solving the task of improving the oil system.
The procedure for calculating the required pumping of oil through the support and a mathematical model of thermal processes in the oil system are often based on equations of thermal balance and heat transmission for a single-phase environment. The models that take into consideration the multiphase nature of the flow hold only for individual geometries. Given the need to blow the seals, the flow in the oil cavity of the GTE rotor support is always two-phase, with a high volumetric share of air. At the same time, due to the rotation of the flow and gravity, there is significant heterogeneity in the distribution of oil and, as a result, in the heat transfer coefficient along the wall of the chamber that forms the oil cavity. Experimental data do reflect these features,

Introduction
The development of aircraft gas turbine engines (GTE) is inextricably linked to increasing the degree of rise in the pressure and temperature of the gas at the inlet to the turbine while reducing the mass and dimensions of the structure, the cost of auxiliary systems. Improving the efficiency of GTE and developing new energy-efficient technologies in the engine industry inevitably sets the tasks of improving the oil systems, which is implemented by reducing the supply of oil to friction nodes in the supports of the GTE rotors and by the rational design of the support nodes. Insufficient supply leads to the intensification of wear and overheating of support elements, destruction of bearings.
A promising direction in this regard is to model an oil system at the stage of its design, where the thermal-hydraulic processes in the supports of GTE rotors are the

CFD MODELING OF MULTIPHASE FLOWS IN THE GAS TURBINE ENGINES OIL CAVITIES
in a stationary statement for incompressible phases. The computations involved the Realizable k-ε turbulence model and the implicit time-dependent sampling.
Calculating the trajectory of oil particle movement, taking into consideration the use of steady and unsteady solvers with different time steps, was reported in [12]. The authors of work [13] used the VOF method and the SST k-ω turbulence model with Turbulence Damping, as well as the Compressive Interface Reconstruction calculation scheme with the explicit second time sampling scheme, implemented in the ANSYS Fluent, to model processes in the bearing chamber. The application of such a scheme significantly accelerated calculations while maintaining the required accuracy. A similar set of models and settings, but applied to the modeling of oil leaks through the labyrinth seals of the bearing chamber, was proposed in [14].
Paper [15] reports an approach to simulate film behavior at the bearing chamber walls, combining the VOF model and the Euler thin-film model (ETFM). Such a model is used for regions with a particularly thin film (less than the height of the mesh cell), and for other regions, where the mesh's resolution is sufficient to describe the thickness of the film, a standard VOF model is used. The two models are interconnected by the source terms, which are part of the mass and momentum equations, and the criterion that determines the applicability of a particular model. This criterion was the volume of a liquid phase in the cell.
Thus, to date, there is no universal procedure for describing multiphase currents in the oil cavities of the GTE rotor supports. Each separate problem implies individual approaches to the solution, and the construction of a mathematical model implies the introduction of many assumptions. The issue of choosing the structure of a mathematical model, taking into consideration the capabilities of modern CFD packages, remains open, which is why research in this direction is still relevant.

The aim and objectives of the study
The aim of this study is to substantiate the procedure for constructing and improving a mathematical model of the thermohydraulic processes in the oil cavities of the GTE rotor supports. This would make it possible to improve the accuracy of the calculation of the temperature state of the elements of the GTE rotor support, to numerically investigate the applicability of those known and to derive the new correlation dependences for the mean value of a heat transfer coefficient. A possibility to numerically investigate the impact of the geometry, the rotor rotation frequency, and the flow rate of phases for heat output in the support's oil cavity could reduce the volume of experimental-finishing operations.
To accomplish the aim, the following tasks have been set: -to justify the choice of the structure of an oil-air flow model in the oil cavity of the GTE rotor support based on the analysis of the workflow; to build a three-dimensional CFD model of thermal-hydraulic processes in the oil cavity using the information on flow distribution and heat exchange given in the scientific literature; to analyze the impact of common approaches, models of individual processes, and boundary conditions, on the results of calculations and their alignment with experimental data. but studies have been performed for oil cavities in the supports of simplified geometry. The generalizing dependences are rare and their applicability for oil cavities of the rotor supports of an aircraft GTE needs to be checked given significant differences in the geometry and mode parameters.
The heat transfer coefficient is a key parameter when calculating the temperature state of the rotor support elements. The relationship between the thermal and hydrodynamic processes, the complexity of the chamber geometry that forms the oil cavity and the boundary conditions require that its determination should employ simulation techniques based on the mechanics of multiphase environments and computational fluid dynamics (CFD). Research conducted by Safran, Rolls-Royce, and others companies testifies to the relevance of this task and such an approach to modeling.

Literature review and problem statement
The use of computational fluid dynamics (CFD) methods has expanded the possibilities of research into flow distribution and heat exchange in oil cavities. Since 2000, the University of Nottingham Technology Centre has reported many results of the numerical simulation of an air and oil flow in the aircraft engine bearing chamber. The application of the commercially available CFD CFX software suite for the numerical study of turbulent airflow and the trajectories of oil droplets in the experimental chamber was described in [1]. The Navier-Stokes's time-averaged equations and the k-ε turbulence model were used to describe the movement of airflow. The movement of droplets was calculated using a Lagrange approach. A similar approach was implemented in [2,3] when studying a two-phase air-oil flow in the bearing chamber of the Rolls-Royce Trent aircraft engine.
The combined approach by Euler and Lagrange to calculate the flow-distribution in the bearing chamber was proposed in works [4,5]. The calculations were carried out in the ANSYS Fluent software package. The two-phase flow of air and oil droplets was described using then Euler's approach, implemented in the VOF (Volume of Fluid) model of the homogenous multiphase current [6]. All phases were considered as solid, regardless of their actual morphology. The authors introduced the concept of a volumetric fraction of the phase as another additional flow parameter. For all phases, one system of equations is solved, including the equations of maintaining momentum and the transfer of the volumetric share of each phase.
Paper [7] describes a one-dimensional hydrodynamic model that makes it possible to determine the average velocity of airflow in the chamber, the tangential stress at the surface of the film, and the distribution of the oil film along the circumference of the bearing model chamber. It combines the models of film motion [8,9] and airflow motion in the chamber's circular space [10]. These models are interconnected by the shear stress from the air side, used in the film motion model, and the velocity of film motion at the air-film contact boundary used to close the interrelations of the airflow motion model.
The results of a numerical simulation of the two-phase flow of liquid-gas in a cavity with a rotating shaft are reported in [11]. The calculations were carried out in the ANSYS Fluent software package using the Euler homogeneous multiphase flow model VOF, taking into consideration the forces of surface tension and mass forces. The problem was solved

Materials and research methods
The object of our study is the support's oil cavity (Fig. 1), described in works [16,17]. Given the limited amount of publicly available information on the results of experimental studies in oil cavities and about the detailed conditions of the experiment, and in order to verify the CFD model that we constructed, chamber II was examined ( Fig. 1), for which the results of the experiments are available [16,17].

Fig. 1. Experimental bearing chamber
The heated oil enters, through a roller bearing, the experimental chamber II, formed by a rotating shaft and a fixed body. Through a labyrinth seal, the air is fed into the chamber; the air has the same initial temperature as the oil. A mixture of oil and air comes out of the chamber through the venting (atop the chamber) and scavenge (bottom of the chamber) pipelines.
When constructing a mathematical model, a turbulent stream consisting of two mutually insoluble phases with the same pressures is considered; in it, the carrier phase is air, the dispersive phase is oil drops. The liquid phase is exposed to two main forces: gravity and viscous friction. For a multiphase heterogeneous flow, conservation equations are applied separately for each phase.
The continuity equation: The momentum equation: The energy equation: The volumetric density of phase i is calculated from 0 .
For equations (1)-(4), the following designations are adopted: λ is the full enthalpy, density, dynamic viscosity, and thermal conductivity of phase i, respectively; U i is the phase i velocity vector; J ji is the mass flow from phase j to phase i in a unit of volume; F Ar is the lift force; F a is the force of the attached masses; S mi , S di , S Ei are the source terms of the phase i mass, momentum, and energy, respectively; characterize the interphase exchange of momentum and the thermal interaction between phases i and j, respectively; α i is the volumetric share of the phase. Their sum The above equations describe patterns in the movement and heat exchange of an oil-air mixture inside the chamber and an oil film, which is formed on the walls. For the case of a separate description of the main stream and the film, additional interaction conditions at their interface boundary are recorded.
The examined bearing chamber and the boundary conditions are schematically shown in Fig. 2.
When modeling a multiphase flow, the boundary conditions that were assigned included pressures at the outlet from the venting and scavenge ports, the flow rates and temperatures of air and oil streams entering the chamber, temperatures at the surface of the stator, bearing, and shaft, the shaft rotation frequency, the diameter of droplets reflected from the bearing, heat flux through the transparent wall.
A geometric model of the chamber (Fig. 3) was built to carry out numerical modeling. In order to reduce computational time costs while maintaining the acceptable accuracy of calculation, we used a structured grid with a minimal step near the walls and a maximum step within the main volume. Based on the results of studying the mesh convergence, as well as the effect exerted on the result by the quality and size of the grid, further research involved an option with a total number of elements of 0.7 million, where the minimum cell size is 0.02 mm and the maximum -7 mm (Fig. 4). The numerical simulation was carried out using the Ansys CFX and Fluent software suits, 2019 R2 version.

1. Justification for choosing the structure of the mathematical model of an oil-air flow
As regards the support's oil cavity, the phase distribution and a heat transfer coefficient from the wall are the main characteristics. To determine them, it is necessary to solve conjugal thermohydraulic problems in a three-dimensional statement for a drip-air flow and the near-wall oil film under the boundary conditions that depend on the geometric and mode parameters.
The mechanics of heterogeneous environments employ both the Euler's continual approach and Lagrange's trajectory approach. When using the latter, the velocity field of a disperse phase is determined without taking into consideration the impact on the carrier phase, the conditions of heat exchange of the flow with the elements of the chamber are not accurately described. As a result, the temperature and velocity of the droplets are not determined correctly enough, which further affects the parameters of the near-wall oil film and, consequently, the accuracy of the calculation of a heat transfer coefficient. Given this, further modeling of the drip-air flow was carried out using the Euler's approach for both phases.
Based on the analysis of the scientific literature [11,15,[18][19][20], preliminary calculations [21][22][23][24], and recommendations on the application of models in Ansys [25], 12 structures of a three-dimensional CFD model were used for this study. The best consistency with experimental data was shown by structures based on the application of the Inhomogeneous (Ansys CFX), Euler and VOF (Ansys Fluent) models.
The Ansys CFX software suite implements the Inhomogeneous model as part of the Particle Model, which is designed to describe the behavior of a two-phase flow containing a continuous and dispersed phase. Based on the Eu-ler's approach, we modeled a heterogeneous monodispersed oil-air flow with the condition of inversion of the structure of the two-phase flow in the near-wall area from drip to bubble, which made it possible to exclude from consideration a mathematical model for phase three, the oil film. The area of the interphase surface in the unit of the mixture volume is determined from the following formula: 6 , At high volumetric share values α j , there may be a violation of the condition of the j phase existence as a dispersal phase, which corresponds to the inversion of the structure of the two-phase flow. In this case, the Particle Model imposes a limit on the volumetric share α max =0.8 in accordance with the algorithm: The condition α min =10 -7 excludes the interruption of counting when the volumetric share of the disperse phase is approaching zero.
The Ansys Fluent used the VOF and Euler models in stationary and non-stationary statements. The VOF model solves three motion equations based on phase velocity averaging, a continuity equation for one phase, and an equation of conservation of a volumetric share. When phase velocities differ significantly, the accuracy of the model is compromised. The Euler model solves six motion equations (three for each phase), a continuity equation (for one phase) and a conservation equation for a volumetric share. The peculiarity of the Euler model is that the motion equations are solved for each phase separately and linked together by the interphase interaction factor.

2. Developing a three-dimensional CFD model of the thermal-hydraulic processes in an oil cavity
To date, there is no unified method for analyzing multiphase flows. Each individual problem requires an individual approach that takes into consideration the specificity of a workflow in the device and the basic parameters that need to be defined.
Based on the analysis of literary data and test calculations, we have chosen a Pressure-Based solver and the Realizable k-ε turbulence model.
Transfer equations for the turbulence kinetic energy k and energy dissipation rate ε are recorded in the following form: Here, μ t is the turbulent viscosity; σ k and σ ε are the turbulent Prandtl numbers; G k and G b is the gain in the turbulence kinetic energy at the expense of a gradient of the averaged velocity and lift force, respectively; Y M is the energy change at the expense of turbulent flow compressibility; S k and S ε are the source terms; C 1 , C 2 , C 3 , are the constants of the model; ν is the kinematic viscosity.
The near-wall functions Scalable Wall Functions and Enhanced Wall Treatment were used to correctly describe a flow near the wall.
The results of modeling the distribution of oil and air flows in the chamber are shown in Fig. 5-7. The drops of oil under the influence of centrifugal forces, and the airflow directly, accumulate on the wall of the chamber, forming a film. According to the results of the calculation (Fig. 5, 7), the thickness of the oil film on the circumference of the chamber is variable. The velocity of film motion is determined by the force of gravity and shear at the film-air interphase.

3. Analysis of the impact of approaches, models, and boundary conditions on the calculation results
Our analysis of literary sources and our calculations have shown that the processes of heat-mass exchange and hydrodynamics in the oil cavities of the rotor support can be described through different models, their settings, and the types of the solver. And given the fact that each subsequent version of the CFD software suite is complemented by new features, the issue of the degree of applicability of a particular model becomes even more acute.
At the preliminary stage of the simulation, those structures and models that performed incorrectly or caused problems with convergence were discarded. The most stable solution was provided by the boundary conditions of "Mass flow inlet" at the oil and air inlets, as well as "Pressure outlet" at the outlets from the chamber.
The results of the comparison of calculations based on the constructed model involving the steady and non-steady solver, as well as two fundamentally different models, Euler and VOF, for a control cross-section of 15° (relative to the venting pipeline in the direction of shaft rotation) [22] of the bearing chamber [16,17] are shown in Fig. 8-10.  9 shows that Taylor's toroidal vortexes form in the central part of the chamber. A similar pattern is observed in the case of the flow motion between two concentrically placed cylinders when the outer cylinder is stationary and the inner cylinder rotates. In this case, the circumferential flow velocity changes its value from zero at the surface of the stationary cylinder to the velocity of rotation of the moving wall.
In order to select the structure of the CFD model that suits best the problems within the examined class, we compared the estimated and experimental data [16,17] Fig. 11.
The results of the comparative analysis show a qualitatively similar pattern of the change in a heat transfer coefficient, which indicates the proper choice of the meaningful model when calculating the flow distribution and heat exchange in the bearing chamber. In the range of 0...10 mm lengthwise the chamber in the examined cross-section, the Euler model involving a non-steady solver (the Fluent calculation) more accurately describes the flow behavior and a heat transfer coefficient. Over this region, the margin of error does not exceed 10 %. The maximum error of calculations does not exceed 30 %. Electronic copy available at: https://ssrn.com/abstract=3708308 The impact of gravity (Fig. 12) was also investigated as part of our work. The lowest heat transfer coefficient values when gravity force is considered are observed in the cross-sections with angular coordinates 195-240° relative to the position of the venting pipeline. Increasing the thickness of the film in this region increases the thermal resistance of the boundary layer and, consequently, the heat transfer coefficient drops. In the range of angular coordinates of 270-315°, where the formation of the film is difficult and its thickness is insignificant, the heat transfer coefficient reaches a maximum value. The chart in Fig. 12 shows that disregarding the force of gravity can lead to a difference in the heat transfer coefficient of up to 50 %.

Discussion of the results of studying the thermalhydraulic processes in the oil cavity of the GTE rotor support
Our results of numerical modeling show that the proposed model properly describes the thermal-hydraulic processes in the oil cavity of the GTE rotor support. It should be noted that the accuracy of the calculation based on the proposed model can be improved, in particular, by making the grid smaller and by even the more detailed mathematical notation of the processes in the near-wall area. However, this will inevitably lead to a significant increase in the time of computation, which will require even greater computer resources and will inevitably go beyond engineering calculations. Thus, for a given problem, the average time of computation using the VOF model was about 700 hours (computer with 16 cores), for the Euler modelabout 1,300 hours.
The uncertainties associated with describing the behavior of the oil film by setting the parameters for a drop flow from the bearing and conditions of interphase interaction in the near-wall area make it difficult to use existing models of the workflow in the oil cavity of the GTE rotor support. As the results of CFD modeling (Fig. 5) show, the thermal-hydraulic processes in the right and left parts of the chamber are not symmetrical. This is primarily due to a change in the direction of the shear (shear stress) force at the interphase surface film/air+drops relative to the gravity force.
Due to the impact of gravity, the film moves downwards. At the same time, the shaft in its rotation involves the air that got into the chamber through the labyrinth seals into circular motion (Fig. 6, 7). Where the directions of airflow and the film are the same (Fig. 5, right side), the velocity of the latter increases while the thickness decreases. In the other part of the chamber, where the air flows and film are diverted (Fig. 5, left side), there is braking of the film with an increase in its thickness. In addition, in the left part of the bearing chamber, there is an accumulation of oil near the venting port where it enters the channel. This result is confirmed by experiments [26].
The comparison of the results of CFD calculations (Fig. 8-10) demonstrates, in general, qualitatively close values of the main characteristics of the flow for different models and solvers. However, there are also differences, which can lead to a significant error in determining the heat transfer coefficients from the wall to the two-phase flow in the chamber. Thus, when comparing the temperature fields (Fig. 8) for different calculations, the differences in temperature in some cross-sections reached 8 K.
When comparing velocity vectors (Fig. 9), certain discrepancies in the size of vortex zones and their location were identified. In the cross-sections near the oil inlet, there was a formation of additional local vortexes. At similar points of space for different calculated cases, differences in the velocity of individual phases reached 7 m/s, and mixtures -5 m/s (Fig. 10). The greatest discrepancy was observed in the near-wall regions and places of sharp flow turn, which can be explained both by the complexity of the physics of the processes themselves and by the imperfection of the mathematical apparatus to describe them. The comparison of calculated and experimental values for a heat transfer coefficient (Fig. 11) reveals that the greatest discrepancy is observed at the edge of the wall far from the oil inlet (a distance of 18-20 mm). This can be explained by the presence of vortex and stagnant zones in this part of the chamber, the imperfection of the existing mathematical apparatus to describe such zones, and the inability to fully recreate the conditions of the experiment in the model. Gravity force has a qualitatively important impact on the distribution of the film. Due to its asymmetric distribution (Fig. 5), the complex movement of air (Fig. 6) and oil droplets (Fig. 7), there are different values of thermal resistance of the boundary layer and, hence, the heat transfer coefficient. This fact is confirmed by the results of calculations based on the model that we developed (Fig. 12).
It should be noted that in the proposed model the diameter and velocity of oil droplets flying out of the bearing are set as the boundary conditions. Therefore, the dependence of these values on the regime parameters is not fully taken into account. The issue of setting the boundary conditions on the surface of the bearing and temperature field of the parts adjacent to the study area is also not completely resolved. This is especially true for those oil cavities in which a bearing simultaneously serves as an inlet and outlet for oil fed into Heat transfer coefficient, W/m 2 K θ, degree considering gravity force neglecting gravity force Electronic copy available at: https://ssrn.com/abstract=3708308 the chamber by nozzles. The answers to these questions can be obtained mainly from experimental studies of actual GTE bearing chambers.

Conclusions
1. The choice of the structure of a working process model in the bearing chamber has been substantiated, based on the concept of interpenetrating continuums, containing the algorithms that take into consideration the formation of a near-wall oil film and its interaction with the adjoining drip-air flow.
2. To determine the phase distribution and a heat transfer coefficient in the oil cavity, a three-dimensional CFD model of thermohydraulic processes was built based on the Inhomogeneous and Particle models implemented in the software Ansys CFX, as well as the Euler and VOF models, given in the Ansys Fluent package. The model implies the use of a steady and non-steady Pressure-Based solver.
3. The impact of approaches and models on the results of the calculation has been analyzed. The best consistency with the experimental values of a heat transfer coefficient was shown by the Euler model using a non-steady solver, the Realizable k-ε turbulence model, and the Enhanced Wall Treatment near-wall function. It is shown that the uneven flow distribution in the bearing chamber leads to a variable value of the heat transfer coefficient, which has a significant effect on the heat flux into the oil cavity of the rotor support, on the temperature state of the oil and support elements. Gravity significantly affects the formation of an oil film and must be taken into account in the calculation.