## Abstract

In the case of aero-engine, thin lubricating film servers dual purpose of lubrication and cooling. Prediction of dry patches or lubricant starved region in bearing or bearing chambers are required for safe operation of these components. In this work, thin liquid film flow is numerically investigated using the framework of the Eulerian thin film model (ETFM) for conditions, which exhibit partial wetting phenomenon. This model includes a parameter that requires adjustment to account for the dynamic contact angle. Two different experimental data sets have been used for comparisons against simulations, which cover a wide range of operating conditions including varying the flowrate, inclination angle, contact angle, and liquid–gas surface tension coefficient. A new expression for the model parameter has been proposed and calibrated based on the simulated cases. This is employed to predict film thickness on a bearing chamber which is subjected to a complex multiphase flow. From this study, it is observed that the proposed approach shows good quantitative comparisons of the film thickness of flow down an inclined plate and for the representative bearing chamber. A comparison of model predictions with and without wetting and drying capabilities is also presented on the bearing chamber for shaft speed in the range of 2500 RPM to 10,000 RPM and flowrate in the range of 0.5 liter per minute (LPM) to 2.5 LPM.

## 1 Introduction

In industrial situations where a liquid film is used for thermal management, ensuring an even undisrupted film over a solid substrate is critical for efficiency and for some applications is a requirement for safe operation [14]. When a liquid film is disrupted, it is known as a partially wetted film. Dry patches are formed on the solid substrate resulting in areas with no liquid film or cooling effect. Partial wetting films are observed in the case of aero-engine bearing chambers at low shaft speeds. Lubricating oil, after being shed from bearings, is collected at the walls of the chamber forming a thin film which is driven by gravity and air shear. Depending on the shaft speed and lubricant flowrate, this film may cover the chamber walls partially or fully as discussed by Kurz et al. [5]. Partially wetting films are also observed in the wake of geometric obstructions such as nuts and bolts, Eastwick et al. [6].

Lubrication oil also serves as a coolant. Therefore, there is a need to accurately predict when a liquid film will become disrupted and the format/size of any dry patches to ensure that these conditions can be avoided or systems can be designed to cope with these scenarios. In addition, it is necessary to achieve a good understanding of the thin film flow behaviors under these conditions and to validate the numerical models used in the design and analyses of relevant engineering systems. Particularly, reliable and computationally economic methods will help engineers and designers to predict the performance and to optimize their operational parameters in cost-effective time.

A lot of work has been performed in the past both experimentally and numerically to understand thin film flows over a solid substrate. A falling film over an inclined flat plate can show complex flow behaviors and exhibit different flow regimes and surface wave instabilities, Nusselt [7]. The stability of thin film flow and formation of surface waves have received great attention demonstrated by a large amount of literature work. For example, recently, Denner et al. [8] studied the hydrodynamics of solitary waves on falling liquid films using the volume of fluid (VOF) and Eulerian thin film model (ETFM) methods and compared the results against experimental measurements. The VOF model compared very well with experimental measurements and the ETFM results showed accurate prediction of the film height even for high Reynolds number conditions. This shows the ability of this approach to capture complicated thin film flow behaviors.

Other complex thin film flow phenomena are related to partial wetting and transition from the film flow regime to rivulet flow or droplets depending on parameters such as flowrate, inclination angle, and surface finishing. Huppert [9] was the first to study such instabilities of viscous fluids flowing down an inclined plate, noting the break-up of flow into fingers. Following that work, other experimental studies were performed to investigate different flow pattern formation for fluids with different wetting properties and to characterize the instabilities that develop at the contact line [10,11].

Troian et al. [12] theoretically studied the instability of a liquid film flowing down an inclined plate and found that the instability is the largest at the advancing front, which shows capillary ripples near the contact line. More recently, several studies have been conducted to investigate the film breakup and rivulet flows both experimentally and numerically with a focus on the VOF method to predict the wetting flow behavior. Lan et al. [13] studied a partially wetting falling film flowing down an inclined plate both experimentally and numerically using the VOF method and with varying the flowrate, surface tension, contact angle, and inclination angle. The film thickness and widths were measured using a laser focus displacement instrument. The comparisons of numerical simulation and measurements showed good agreement. The effect of liquid film flow rate and surface treatments were numerically studied in Iso-and Chen [14] using the VOF framework. The comparisons with existing experiments and validations showed that the method is capable of predicting the transition between film flow and rivulets. Bonart et al. [15] performed experimental investigations on the effect of surface tension, viscosity, inclination angle, and mass flow on the interfacial area of rivulets and developed a correlation based on Reynolds and Kapitza numbers.

Three-dimensional VOF simulations were conducted by Singh et al. [16] to calculate the interfacial area of rivulets for a wide range of varying parameters such as inlet size, inclination angle, contact angle, flowrate, and solvent properties. A scaling proposal for the interfacial area as a function of Kapitza number and contact angle was introduced. In another work, Singh et al. [17] studied rivulets stabilities and breakup on an inclined plate using VOF. The CFD simulations compared well with experiments for liquids with various viscosities, and they developed a criterion that determines rivulet and droplet flow regimes based on Weber number as a function of Kapitza number. From the above, it can be concluded that VOF simulations can generally offer sufficiently accurate results for predicting partial wetting phenomena. However, it is important to consider the CPU time needed to perform these simulations, which is prohibitive for aero-engine bearing chambers.

Meredith et al. [18] studied a partially wetting film with rivulets both experimentally and numerical simulations using thin film model assumptions. The contact angle shear stress was derived based on Young's law and implemented within the ETFM in the open-source openfoam CFD solver [19]. The model was able to predict the transition from film flow to rivulets; however, the inlet film mass flow rate in the simulations was modified to match the flow regimes in the experiments. Calibrating a numerical model by altering inlet mass flow rate is not tenable for an efficient predictive tool. In later work by Martin et al. [20], an ETFM was implemented in OpenFOAM, and it was used to study different types of gravity-driven thin film flows. The simulation results showed good agreement for wavy film flow in terms of film thickness, wave speed, and their energy spectra. Very limited qualitative comparisons of rivulet flow for two selected cases, from the experimental work of Silvi and Dussan [10], were presented using a prescribed advancing contact angle. Singh et al. [21,22] proposed coupling of ETFM with VOF model to simulate a bearing chamber. The coupled approach was successfully implemented in these studies, but the quantitative comparison of predicted results was not in good agreement with the experimental results. These studies also highlighted the scope of further improvements in the ETFM.

Within the open literature, there is a lack of comprehensive and quantitative validations, and assessment of the ETFM and is extended for different operating conditions, which motivates this work. The thin film model (ETFM) is of particular interest because it offers a computationally cheaper alternative to surface tracking methods for thin film flows. Typically, the film thickness in aero-engine bearing chamber oil films is less than 1mm, Eastwick et al. [6]. Highly refined meshes are needed for surface tracking techniques to capture the sharp liquid–gas interface and its interaction with the gas flow in three-dimensional simulations. In addition, the time stepping requirement for such fine meshes is prohibitive as reported by Bristot [23] and Kakimpa et al. [24]. However, since the depth-averaged equations are solved in the ETFM, there is no requirement for a fine mesh at the liquid–gas interface. Therefore, if improvements can be made in ETFM predictions, significant cost saving can be made.

In this work, a liquid thin film flow is numerically investigated using the framework of the ETFM, and this is extended to include a new term in the model with parameters calibrated using an extensive dataset. A partial wetting model has been implemented and integrated with ETFM in ansysfluent within the framework of the Clean Sky 2 AERIS project, involving the University of Nottingham, ANSYS (ANSYS Inc., Sheffield, UK), and Rolls-Royce (Rolls-Royce plc, Derby, UK). The project is developing multiphase flow models for simulations of air-oil flow in aero-engine bearing chambers, where partially wetted film flow is relevant at low shaft speeds. The developed model is available in ansysfluent v19.2 [25], which is used in this work. An assessment of the model's capability to predict the partial wetting phenomena for film flows and comprehensive validations against experimental data are presented. The optimized model is employed to predict film thickness on a bearing chamber for shaft speeds in the range of 2500 RPM to 10,000 RPM and liquid flow rates in the range of 0.5 liter per minute (LPM) to 2.5 LPM.

## 2 Problem Description

The bearing chamber considered in this study is adopted from the experimental rig of Chandra et al. [26]. This consists of a chamber, film-generator, shaft, and sump as depicted in Fig. 1(a). The nondimensional sizing of the chamber is demonstrated in Fig. 1(b). Water is supplied to the inner surface of chamber via a film-generator. Water enters from the lip of the film-generator onto the inner surface of chamber at 37 deg clockwise from the top dead center (TDC) of the chamber and flushed to the sump 156 deg clockwise from TDC. The computational model used in this study mimics all the dimensions of the experimental chamber. A film mass-flux is defined at the lip of film generator corresponding to the experimental the flow rate.

Fig. 1
Fig. 1
Close modal

### 2.1 Mathematical Model.

In the derivation of the ETFM, the wall normal velocity is assumed to be zero and in addition, convection takes place in the wall tangential directions only. The assumption of a thin film is valid when the thickness of the film with respect to the transversal directions and with respect to the radius of curvature of the surface is much less than 1. In this way, it can be assumed that the film flow is parallel to the wall, and film properties do not vary across the film thickness. Thus, the depth-averaged continuity equation for a Newtonian two-dimensional (transversal direction) incompressible film flow over a wall with a free surface exposed to a gas is given by
$∂ρlh∂t+∇s·(ρlhul)=Sm$
(1)

where ρl is the liquid density, h is the film thickness, $∇s$ is the surface gradient operator, ul is the velocity of the liquid, and Sm is the mass source term due to, e.g., droplet interaction with the film or other surface exchange phenomena. This study is conducted at low liquid flowrate and isothermal conditions. Hence, mass source term is assumed to be zero for all the investigated cases.

The momentum conservation of thin liquid film flow is given by
$∂ρlhul∂t+∇s·(ρlhulul+DV)=−h∇sPL+ρlhgs+32τfs−3μlhul+τθw+Si,mom$
(2)
where the transient and inertia effects are taken into account by the first and second terms on the left-hand side, respectively. A correction tensor, DV, is added to the inertia term in the parentheses to account for the effect of nonuniform velocity distribution across the film thickness. In addition, it was found that this correction provides solution stability for shock and pool flow regimes for thin film rimming flow [24,27]. A quadratic velocity profile assumption is more consistent for the condition of a laminar flow of a liquid film over a solid wall. Hence, in this study, a quadratic velocity profile is assumed. The variable PL on the right-hand side of Eq. (2) is defined in Eq. (3)
$PL=Pg−ρh(n·g)−σ∇s·(∇sh)$
(3)

which accounts for the gas pressure, gravity force normal to the wall, and the surface tension (σ) force based on the film thickness curvature, respectively. The gravity force in the wall tangential direction is represented by the second term on the right-hand side in Eq. (2). Shear forces on the gas-film $(τfs)$ and film-wall interfaces are given by the third and fourth terms on the right-hand side in Eq. (2). For the cases simulated in this study, the gas velocity is very small, and the gravity force in the tangential direction is the main force that drives the flow. So, the gas-film shear force is considered of negligible effect in this study. No momentum sources, Si,mom, are employed in this study.

The contact angle shear stress applied at the contact line between dry and wet areas has the effect of restricting the liquid film from expanding over the solid substrate. This enhanced model is implemented in ansysfluent 19.2 [25]. Following the derivation in Meredith et al. [18], the contact angle shear stress is given by
$τθw=βσ(1−cosθm)∇sw$
(4)
where w is the film wet area fraction and takes a value of zero in dry regions and a value of 1 in wet regions,$θm$ is the mean contact angle. The contact angle is the angle between a liquid–vapor interface and a solid surface. The surface is considered dry if the wall film thickness is below a critical value, hc, which is chosen to be 1.0 × 10−10 m. The implementation ensures that this shear stress is applied only at the contact line. The parameter, ß, is an empirical model parameter that needs optimization to reproduce the experimental behavior. Either a constant static contact angle can be provided as the input or a dynamic contact angle can be assumed. It should be mentioned that the contact angle shear stress given by Eq. (4) is derived from a force balance for a single static droplet on a horizontal surface, Meredith et al. [18]. In real flow conditions, the contact angle can differ from the static one [28,29], and in order to account for this effect in the numerical model, the contact angle can be allowed to vary randomly following a Gaussian distribution, given by
$f(θ)=12πδ2e−(θ−θm)2δ2$
(5)

where δ is the relative standard deviation (SD) with respect to the mean contact angle (θm).

This approach of simulating dynamic contact angle might not be sufficient in wetting and drying phenomena under flow conditions when the contact line is in motion and advances over a substrate surface with different inclination angles. To account for this effect, an additional model empirical parameter, ß, is included in Eq. (4) and requires adjustments for each individual case to reproduce the observed behavior in the experiments. As will be demonstrated, the flows are particularly sensitive to this parameter, and this paper seeks to address this issue through proposing a new correlation which is calibrated on an extensive dataset.

### 2.2 Validation Test Cases.

The experimental data of two different studies are used to validate the numerical methodology, and thereafter, model parameters are optimized. The geometrical setup for the phenomena under investigation is depicted in Fig. 2 following Meredith et al. [18] and Lan et al. [13]. In the numerical model, a thin film originates from an inlet boundary at the top and flowing down under the effect of gravity over a plane surface with an inclination angle ϕ. The inlet boundary has a width smaller than the width of the channel, allowing for the development of a contact line at the boundaries or rivulets in the case of low mass flowrate conditions. The first reference case (computational domain (CD)1) considered in this work is taken from the experimental campaigns performed by Meredith et al. [18], which provide data for a falling water film with varying flow rate and inclination angle exhibiting partial wetting. The second reference experimental case (CD2) considered for validation is taken from Lan et al. [13] where the behavior of a gravity-driven thin liquid film was studied experimentally and numerically, varying the flow rate, inclination angle and surface tension of the working fluid. Lan [13] also provides detailed information about the film thickness and its lateral distribution over the substrate surface. The dimensions of the computational domain marked in Fig. 2 for both test cases are given in Table 1.

Fig. 2
Fig. 2
Close modal
Table 1

Dimensions of the adopted computational domain for different cases in mm

Test caseLWHwshs
CD1122061010051050
CD2550101.620.3276.26.35
Test caseLWHwshs
CD1122061010051050
CD2550101.620.3276.26.35

### 2.3 Operating Conditions.

In total, eight test cases were simulated on the first CD1 following studies of Meredith et al. [18] and eleven test cases on the second CD2 following Lan et al. [13]. Using these, the new ETFM model term and correlation are proposed. The operating conditions, for the cases reported in the work by Meredith et al. [18] (case 1–8) and Lan et al. [13] (case 9–19), are summarized in Table 2.

Table 2

Summary of operating conditions for the Meredith et al. [18] (cases 1–8) and Lan et al. [13] (cases 9–19) test cases

CaseQ × 105 (m3/s)ϕθmσ (N/m)Re
11.695750.06931.2
25.565750.069109
38.955750.069175
415.905750.069311
53.7690750.06973.6
66.4390750.069126
710.9190750.069214
825.5390750.069500
92.0260400.042266
102.0260300.042266
112.0260170.026266
120.6660300.04286.9
131.1360300.042258
140.6660170.02686.9
151.1360170.026148
162.0230300.042266
172.0290300.042266
182.0230170.026266
192.0290170.026266
CaseQ × 105 (m3/s)ϕθmσ (N/m)Re
11.695750.06931.2
25.565750.069109
38.955750.069175
415.905750.069311
53.7690750.06973.6
66.4390750.069126
710.9190750.069214
825.5390750.069500
92.0260400.042266
102.0260300.042266
112.0260170.026266
120.6660300.04286.9
131.1360300.042258
140.6660170.02686.9
151.1360170.026148
162.0230300.042266
172.0290300.042266
182.0230170.026266
192.0290170.026266
In this table, Q is the volumetric flow rate, ϕ is the inclination angle of the plate measured with respect to the horizontal plane, θm is the mean contact angle, and σ is the surface tension coefficient, which is constant for all the cases. The Reynolds number, Re, is calculated as
$Re=ufhf/vl$
(6)
where uf and hf are the film velocity and film height, respectively, and νl is the kinematic viscosity of the liquid. Approximate mean values for film velocity and thickness can be estimated from Nusselt theory for laminar film flow over a planar surface assuming no wave formation at the interface as described by Nusselt [7]. They read
$uf=ρlgxhf2/3μl$
(7)
$hf=[3νlql/gx]1/3$
(8)

where ρl is the liquid density, μl is the liquid viscosity, νl is the liquid kinematic viscosity, gx is the gravitational acceleration, and ql is the volumetric flowrate per unit width. The low values of the calculated Re number for the studied cases show that the film flow remains laminar for all the considered cases Takamasa and Hazuku [30].

### 2.4 Boundary Conditions and Mesh.

Wall boundary conditions are applied at the lower, side, and upper surfaces (as described in Fig. 2). In the investigated cases, liquid is flowing on the lower wall and hence wall film flow equations are solved on the lower wall only. Wall-film boundary conditions are imposed on the lower wall. Injection of the film mass flow rate is achieved through a mass source term distributed uniformly over the inlet area of the liquid film inlet. This source term is implemented as a user defined function. The film mass is monitored during the simulation, and the results are extracted when a stabilized value of the total film mass is reached.

### 2.5 Numerical Methods and Solution Procedure.

ansysfluent 19.2 [25] is used to solve the ETFM (Eqs. (1)(5)) as a two-dimensional problem along with Navier–Stokes equations for the air phase in three dimensions and isothermal conditions. For the Navier–Stokes equations, the second order upwind scheme is adopted for discretization of the convective terms, and pressure–velocity coupling is achieved by the SIMPLE algorithm [31]. Second order implicit time integration is used for advancing the air flow time, while an explicit differencing scheme is used for solving the ETFM equations. A subtime-step which amounts to 10 equal divisions of the air flow time-step is used in the solution of the film equations. The time-step was adapted to ensure the film Courant number was well below 1. A second order upwind scheme is used to discretize the convective terms of ETFM, and the absolute residual convergence criterion is set to 10−5 for the momentum and continuity equations.

### 2.6 Mesh Sensitivity Study and Numerical Uncertainty.

In this study, a grid sensitivity analysis is carried out for both of the computational domains. The structured meshes are made of hexahedral and orthogonal elements and are created with uniformly distributed grid points along the plate in the transversal directions. For the CD1, three grids, viz., CD1-Mesh1 (406,560 cells), CD1-Mesh2 (698,964 cells), and CD1-Mesh3 (1,205,568 Cells) are investigated. Grids are refined in such a way that the grid refinement factor (r) is greater than 1.3. Similarly, three grids CD2-Mesh1 (77,056 Cells), CD2-Mesh2 (129,724 Cells), and CD2-Mesh3 (231,770 Cells) are investigated for CD2. The film thickness distributions in the transverse direction for the investigated grids are shown in Fig. 3 for CD1 and CD2, respectively.

Fig. 3
Fig. 3
Close modal

The numerical uncertainty arising from discretization is estimated based on the procedure described by Celik and Li [32]. The film thickness distribution in the transverse direction is selected as a parameter for estimating uncertainty. For CD1, the range of local order of accuracy (p) was varied from 0.226 to 5.583 with a global average (pavg.) of 3.046. Likewise, the local order of accuracy for computational domain CD2 was varied from 0.172 to 8.31 with a global average of 3.85. The discretization errors are quantified in terms of grid convergence index (GCI), which is calculated based on the globally averaged order of accuracy (pavg.). The maximum uncertainty in the film thickness because of discretization error was ±17.14 μm for CD1 at z/W =0.12 and ± 10.47 μm for CD2 at z/W =0. Based on the grid dependence study, Mesh2 (CD1-698,964 cells) and Mesh 3 (CD2-231,770 Cells) are selected for both the domains to conduct the parametric study presented in Secs. 2.72.9. The uncertainty in the film thickness arising due to discretization is plotted in Figs. 3(a) and 3(b) for CD1 and CD2, respectively, along with film thickness variation for the baseline meshes.

### 2.7 Validation of the Computational Methodology.

Initially, validation of the cases in Table 2 is presented. Numerical studies were conducted to reproduce all the results of the 19 test cases presented in Table 2 but due to brevity only selected cases are presented here. Comparisons for the liquid film distribution and the flow regimes are shown in Fig. 4 against the experimental measurements by Meredith et al. [18] for inclination angles of ϕ = 5 deg (cases 1, 2) and ϕ = 90 deg (case 8) in Table 2. The parameter values selected in the simulations corresponding to the results in Fig. 4 were ß =1.0 and SD = 10%.

Fig. 4
Fig. 4
Close modal

At low flow rate conditions as in cases 1 and 2, rivulets are predicted in the simulations as in the experiments. It should be emphasized here that the exact shape of rivulets and their paths vary experimentally due to a number of uncontrolled factors such as surface finish variations or prewetting of the solid substrate, and for this reason, the aim is to compare the general trends of the rivulets in the simulations compared to experiments under the same conditions, i.e., the degree of wetting. With the increase of the flowrate and the Reynolds number, the change of the flow regime starts to take place as depicted in Fig. 4(c), where a continuous film is predicted for test case 8.

The application of the model in these three cases shows that different flow regimes can be predicted reasonably well using the same model parameter given as ß =1.0 without the need for model corrections. Thus, the enhanced ETFM can reproduce the experimental behavior with greater accuracy than that reported in Meredith et al. [18], where changes of the input mass flow rates for the simulations were needed to match the flow regimes in the experiments for each case. More quantitative validations are presented in Fig. 5, considering comparisons against the experimental dataset of Lan et al. [13]. A comparison of the simulation results against experiments are given for the test cases 14 and 15 (θm = 17 deg and ϕ = 60 deg). The model empirical parameter, ß =6.0 is adopted based on optimization for both cases. It should be noted that the results of ETFM model are comparable with the VOF results. Thus, the model shows ability to predict the effect of the increase of the inlet flow rate manifested by higher width and thickness of the developed film at much lower computational cost.

Fig. 5
Fig. 5
Close modal

A fair agreement for the comparisons with the experiments is achieved for both the computational domains though small discrepancies are present. Some of these discrepancies can be attributed to experimental effects reported in Meredith et al. [18] including fluctuation of the film thickness at the inlet in the experimental setup and presence of surface normal velocities. In addition, on the modeling side, more precise information on the behavior of the dynamic contact angle could help increase accuracy of simulations. Nevertheless, the enhanced ETFM, as shown in this work, can reproduce the experimental behavior with more accuracy than that reported in Meredith et al. [18].

### 2.8 Sensitivity Study on Model Parameters.

A numerical study is conducted to show the influence of model empirical parameter (ß) on the prediction of partial wetting. Case 5 (Table 2) has been simulated assuming values of ß that range between 1 and 5 as shown in Fig. 6. In Fig. 6(a), when ß =1.0, the contact line advances and progressively wets the solid substrate until it reaches the bottom boundary. In addition, an increased film thickness is developed at the boundary contact lines due to the effect of the contact angle shear stress, which acts to pull the interface away from the wall. Increasing ß values to 2.0 (Fig. 6(b), film breakup is predicted at a midway distance between the inlet and the outlet, and rivulets are formed afterward. In this case, the contact line shows high curvatures at the early stages following the flow release and develops into rivulets in a similar way as observed in the early experimental work of Huppert [9]. With further increases of the model parameter, ß, in Figs. 6(c) and 6(d), film breakup takes place at earlier locations, and the predicted rivulets become thinner. In addition, it is observed that the distance between rivulets increases with increasing ß values. From these results, it can be concluded that the contact angle shear stress can have a dominant effect in determining different flow regimes and also in the predicted film thickness. Therefore, calibrating this parameter ß is highly important in correctly predicting the film flow characteristics under different operating conditions using the ETFM approximations. Currently, no guidelines or correlations exist for the setting of this parameter.

Fig. 6
Fig. 6
Close modal

In order to further evaluate the model for predicting the wetting and drying effects of falling liquid films, the experimental data from Lan et al. [13] are used for comparison and assessment of the model and to develop a better understanding about the behavior of the model parameters and the influential variables on the wetting and drying phenomena. These data from Lan et al. [13] offer the possibility of local quantitative comparisons of both film width and thickness at one downstream location.

A comparison of the film thickness and width at a downstream distance of 138.1 mm is illustrated in Fig. 7 for test case 9. The SD is set to zero to match the conditions used in the simulations by Lan et al. [13]. It is shown in Fig. 7 that by adopting certain values for the model empirical parameter, ß, that results of the same quality as VOF [13] can be achieved for both the film thickness and width; however, with much lower computational expense which show the potential of the ETFM to predict the relevant phenomena with sufficiently good accuracy. However, the value for ß has to be established in order to achieve this match.

Fig. 7
Fig. 7
Close modal

Therefore, several sensitivity studies have been performed on the model empirical parameter, ß, for each of the cases shown in Table 2 in order to achieve good agreement with the experimental data, regarding flow regime, film thickness and width. These data were then used to develop a correlation for the optimum value of model empirical parameter (ß) for different conditions. This correlation is discussed in Sec. 2.9.

### 2.9 New Correlation for Model Empirical Parameter.

In order to take into account, the contact angle effect in the ETFM, the contact angle shear stress given by Eq. (4) is added to the momentum transport equation as a source term. This equation was obtained by applying Young's law for a stagnant single droplet on an ideal solid horizontal plane, Meredith et al. [18]. However, in actual flow conditions, the contact angle can deviate significantly from the measured one at equilibrium conditions [28,29]. In particular, a droplet on an inclined plane would take an asymmetric shape with the advancing contact angle larger than the receding one. This suggests that the contact angle force calculated based on the static contact angle would greatly underestimate the effect of the contact angle shear stress for the advancing contact line over an inclined plate. It was shown experimentally in a number of studies [10,11] that the contact angle has a major influence in determining the flow regime and the wetting characteristics of a falling film over inclined plates. This justifies the need to adopt an empirical parameter in the simplified submodel, given by Eq. (4), to account for the effect of the contact line being in motion over an inclined plate. This is particularly important in absence of information about dynamic contact angle behavior for the designated conditions.

Two parameters were identified to play a major role to reproduce the behavior observed in the experimental data; the static contact angle, θm and the inclination angle, ϕ. Correlating optimized ß values for all the simulated cases with these two parameters gave the following new proposed mathematical correlation:
$β=3.69×(sin(ϕ)θm)0.486$
(9)

where θm is calculated in radians. It is found that this expression fits well with all the simulated cases and giving a maximum relative error of less than 10% compared with the optimized value of the empirical parameter. Table 3 compares the optimum ß with the predicted ß. It should be mentioned that the proposed expression is obtained from two entirely different sets of experiments carried out by different groups, yet the model empirical parameter provides good predictions for nearly all the simulated cases.

Table 3

A comparison of optimum ß with the predicted ß for the operating conditions of Meredith et al. [18] (cases 1–8) and Lan et al. [13] (cases 9–19) test cases

CaseOptimum ßPredicted ß% error
110.987−1.239
210.987−1.239
310.987−1.239
410.987−1.239
533.2347.817
633.2347.817
733.2347.817
833.2347.817
944.0942.355
1054.708−5.822
1166.2063.442
1254.708−5.822
1354.708−5.822
1466.2063.442
1566.2063.442
1643.605−9.869
1755.0490.999
1854.751−4.963
1976.656−4.912
CaseOptimum ßPredicted ß% error
110.987−1.239
210.987−1.239
310.987−1.239
410.987−1.239
533.2347.817
633.2347.817
733.2347.817
833.2347.817
944.0942.355
1054.708−5.822
1166.2063.442
1254.708−5.822
1354.708−5.822
1466.2063.442
1566.2063.442
1643.605−9.869
1755.0490.999
1854.751−4.963
1976.656−4.912

## 3 Implementation of Model to Predict Film Thickness on a Bearing Chamber

Bearing chambers of an aero-engine are an example of a complex two phase flow where sealing air interacts with the lubricating oil. Lubricating oil after dispersion from the bearing is collected at the walls of the chamber. A thin film is formed on the inner walls which is driven by gravity and air shear. In order to demonstrate the capability of the numerical model to predict such a complex flow the developed approach is used to predict film thickness on an aero-engine bearing chamber. The bearing chamber considered in this study is adopted from the experimental rig of Chandra et al. [26]. In this experimental study, water is selected as a working fluid since physical properties of water at room temperature are similar to those of Mobil Jet Oil II at engine representative temperatures. Following the experimental work, water is used as a working fluid for the computational study. The thermo-physical properties of water used in this study are: ρ = 998.2 kg/m3, σ = 0.072 N/m, μ = 0.001003 Pa s and θ = 70 deg. The computational methodology is identical to the falling film test cases and not repeated for the bearing chamber for brevity. In this study, averaged inclination angle from the inlet to exit is considered along the circumference in order to calculate optimum ß.

### 3.1 Validation of Numerical Methodology on a Fully Wet Bearing Chamber.

Owing to a complex flow physics inside an aero-engine bearing chamber, the numerical model has to be validated against the experimental results. All the previous validation studies included only a gravity driven film flow. However, in the case of the bearing chamber, the air shear force has a strong influence on the dispersion of the oil film. Hence, the validation studies on the current bearing chamber will show the reliability of the developed approach. Due to the unavailability of the experimental results in partial wetting regime for an aero-engine bearing chamber, it was decided to use the bearing chamber of Chandra et al. [26], which operated in a fully wet film regime.

In the experimental study of Chandra et al. [26], the film thickness on the chamber wall is measured using a laser confocal displacement meter at 75 deg, 105 deg, and 135 deg clockwise from the TDC. The measurement locations are marked in Fig. 8 along with the mesh used in this study.

Fig. 8
Fig. 8
Close modal

The numerical simulations are carried out at a flowrate of 30 LPM and a shaft speed of 15,000 RPM. A comparison of present numerical results with the experimental results of Chandra et al. [32] is shown in Fig. 9. It can be observed that the present numerical results match well with the experimental results. The maximum deviation of 25 μm in the film thickness measurement is observed at 135 deg for a shaft speed of 15,000 RPM. From the presented results, it can be concluded that the numerical model employed is capable of predicting film thickness on the fully wetted bearing chamber with fairly good accuracy.

Fig. 9
Fig. 9
Close modal

### 3.2 Prediction of Wetting and Drying Capabilities of the Developed Model on Bearing Chamber.

Partial wetting of bearing chambers is observed at lower shaft speeds and low flow rates as reported by Kurz et al. [33]. Hence, the present numerical model is tested by conducting a numerical study for three flow rates in the range of 0.5 LPM to 2.5 LPM at a shaft speed of 5000 RPM. In order to highlight the differences in the prediction of film thickness using the standard ETFM and the enhanced ETFM with the wetting and drying capability, both models are employed. Contours of film thickness on the chamber walls are shown in Fig. 10 for both the ETFM and the enhanced ETFM.

Fig. 10
Fig. 10
Close modal

The film thickness predictions of the standard ETFM show that the film thickness increases with the flowrate as depicted in Figs. 10(a)10(c). However, the drying out region is not captured even for the lowest flowrate investigated in this study using ETFM. On the other hand, the enhanced ETFM predicts dry out of the chamber for flow rates of 0.5 LPM and 1.5 LPM, as shown in Figs. 10(d)10(e). A further increase in flow rate yields identical results for both the standard ETFM and the enhanced ETFM as depicted in Figs. 10(c) and 10(f).

The influence of shaft speed on the film thickness is also analyzed for a fixed flowrate of 0.5 LPM. The shaft speed is varied from 2500 RPM to 10,000 RPM. The contours of film thickness on the chamber walls are shown in Fig. 11. It can be observed from this figure that the coverage of the liquid film on chamber wall at higher shaft speeds increases compared to lower shaft speeds. This also clearly shows the influence of air shear stress on the spreading of the liquid film onto the chamber walls and an increase in the number of rivulets.

Fig. 11
Fig. 11
Close modal

## 4 Conclusions

Computational fluid dynamics is used to simulate the partial wetting behavior of falling films adopting the Eulerian thin film model. In order to account for the partial wetting effect, an additional contact angle shear stress is introduced in the ETFM. Two different sets of experimental measurements are used to validate and to assess the model predictions at various conditions. A mathematical expression based on the inclination angle and the contact angle has been derived by matching the simulation results with the experimental measurements. The proposed expression for the modified model matches to optimized values of the model corrections with errors less than 10% for all the considered cases. The new approach has been applied to different cases varying the inclination angle, surface tension, contact angle, and the flowrate. The prediction showed that the model is able to respond correctly to the change of various parameters in terms of film height and width at the measurement locations and also to predict the transition from film flow to rivulet flow at different conditions of flowrate and inclination angles.

The optimized model is employed to perform film thickness simulations in a bearing chamber and found to reproduce film thickness comparable to experimental measurements. The enhanced ETFM model predicted dry out of the bearing chamber for lower flow rates (0.5 LPM and 1.5 LPM) for the investigated range of shaft speed whereas standard ETFM predicted almost uniform film on the bearing chamber for the same range of shaft speeds.

Thus, it can be concluded that the new correlation and modified thin film model show the capability to reproduce the cases investigated in this work and show highly promising results but at a much lower computational cost compared with volume of fluid methods.

## Acknowledgment

The research leading to these results has received funding from the Clean Sky 2 Joint Undertaking under the European Union's Horizon 2020 research and innovation program under Grant Agreement No 724625.

The calculations were performed using the University of Nottingham High Performance Computing Facility and Athena at HPC Midlands+, which was funded by the EPSRC on Grant No. EP/P020232/1.

## Funding Data

• Engineering and Physical Sciences Research Council (Grant No. EP/P020232/1; Funder ID: 10.13039/501100000266).

• Horizon 2020 (Grant No. 724625; Funder ID: 10.13039/501100007601).

## Nomenclature

### Symbols

Symbols

• D =

shaft diameter, m

•
• e =

error, $ea21=|vr1−vr2vr1|$

•
• g =

grid

•
• GCI =

grid convergence index, $1.25ea21r21p−1$

•
• h =

film thickness, m

•
• H =

domain height, m

•
• L =

domain length, m

•
• p =

apparent order of accuracy, $1lnr21(|ln|ε32ε21|+q(p)|)$

•
• Q =

flow rate, m3/s

•
• q(p) =

$ln(r21p−sar32p−sa)$

•
• r21 =

g2/g1

•
• Re =

Reynolds number, $ufhf/vl$

•
• sa =

$1.sign(ε32ε21)$

•
• u =

film Velocity, m/s

•
• w =

film wet area fraction

•
• W =

domain width, m

### Greek Symbols

Greek Symbols

• β =

empirical parameter to account for contact angle force

•
• δ =

standard deviation

•
• ε21 =

h2–h1

•
• θ =

contact angle, deg

•
• μ =

dynamic viscosity, N-s/m2

•
• ν =

kinematic viscosity, m2/s

•
• ρ =

density, kg/m3

•
• σ =

surface tension, N/m

•
• τ =

shear stress, N/m2

•
• ϕ =

inclination angle, deg

Subscripts

• c =

critical

•
• f =

film

•
• l =

liquid

•
• m =

mean

•
• s =

starting

### Acronyms

Acronyms

• CD =

computational domain

•
• ETFM =

Eulerian thin/wall film model

•
• LPM =

liter per minute

•
• VOF =

volume of fluid

## References

1.
Wang
,
T.
,
Faria
,
D.
,
Stevens
,
L. J.
,
Tan
,
J. S. C.
,
Davidson
,
J. F.
, and
Wilson
,
D. I.
,
2013
, “
Flow Patterns and Draining Films Created by Horizontal and Inclined Coherent Water Jets Impinging on Vertical Walls
,”
Chem. Eng. Sci
,
102
, pp.
585
601
.10.1016/j.ces.2013.08.054
2.
Yu
,
Y. Q.
, and
Cheng
,
X.
,
2014
, “
Experimental Study of Water Film Flow on Large Vertical and Inclined Flat Plate
,”
Prog. Nucl. Energy
,
77
, pp.
176
186
.10.1016/j.pnucene.2014.07.001
3.
Eichinger
,
S.
,
Storch
,
T.
,
Grab
,
T.
,
Fieback
,
T.
, and
Gross
,
U.
,
2018
, “
Investigations of the Spreading of Falling Liquid Films in Inclined Tubes
,”
Int. J. Heat Mass Transfer
,
119
, pp.
586
600
.10.1016/j.ijheatmasstransfer.2017.11.142
4.
Giannetti
,
N.
,
Rocchetti
,
A.
,
Yamaguchi
,
S.
, and
Saito
,
K.
,
2018
, “
Heat and Mass Transfer Coefficients of Falling-Film Absorption on a Partially Wetted Horizontal Tube
,”
Int. J. Therm. Sci.
,
126
, pp.
56
66
.10.1016/j.ijthermalsci.2017.12.020
5.
Kurz
,
W.
,
Dullenkopf
,
K.
, and
Bauer
,
H.-J.
,
2013
, “
Influences on the Oil Split Between the Offtakes of an Aero-Engine Bearing Chamber
,”
ASME
Paper No. GT2012-69412. 10.1115/GT2012-69412
6.
Eastwick
,
C.
,
Huebner
,
K.
,
Azzopardi
,
B.
,
Simmons
,
K.
,
Young
,
C.
, and
Morrison
,
R.
,
2005
, “
Film Flow Around Bearing Chamber Support Structures
,”
ASME
Paper No. GT2005-68905. 10.1115/GT2005-68905
7.
Nusselt
,
W.
,
1916
, “
Die Oberflachenkondensation Des Wasserdamphes
,”
Z. Vereines Dtsch. Ing.
,
60
, pp.
541
546
.
8.
Denner
,
F.
,
Charogiannis
,
A.
,
,
M.
,
Markides
,
C. N.
,
Van Wachem
,
B. G. M.
, and
,
S.
,
2018
, “
Solitary Waves on Falling Liquid Films in the Inertia-Dominated Regime
,”
J. Fluid Mech.
,
837
, pp.
491
519
.10.1017/jfm.2017.867
9.
Huppert
,
H. E.
,
1982
, “
Flow and Instability of a Viscous Current Down a Slope
,”
Nature
,
300
(
5891
), pp.
427
429
.10.1038/300427a0
10.
Silvi
,
N.
, and
Dussan
,
E. B.
, V
,
1985
, “
On the Rewetting of an Inclined Solid Surface by a Liquid
,”
Phys. Fluids
,
28
(
1
), pp.
5
7
.10.1063/1.865410
11.
Jerrett
,
J. M.
, and
De Bruyn
,
J. R.
,
1992
, “
Fingering Instability of a Gravitationally Driven Contact Line
,”
Phys. Fluids A
,
4
(
2
), pp.
234
242
.10.1063/1.858351
12.
Troian
,
S. M.
,
Herbolzheimer
,
E.
,
Safran
,
S. A.
, and
Joanny
,
J. F.
,
1989
, “
Fingering Instabilities of Driven Spreading Films
,”
Epl
,
10
(
1
), pp.
25
30
.10.1209/0295-5075/10/1/005
13.
Lan
,
H.
,
Wegener
,
J. L.
,
Armaly
,
B. F.
, and
Drallmeier
,
J. A.
,
2010
, “
Developing Laminar Gravity-Driven Thin Liquid Film Flow Down an Inclined Plane
,”
ASME J. Fluids Eng.
,
132
(
8
), p.
081301
.10.1115/1.4002109
14.
Iso
,
Y.
, and
Chen
,
X.
,
2011
, “
Flow Transition Behavior of the Wetting Flow Between the Film Flow and Rivulet Flow on an Inclined Wall
,”
ASME J. Fluids Eng
,
133
(
9
), p.
091101
.10.1115/1.4004765
15.
Bonart
,
H.
,
Marek
,
A.
, and
Repke
,
J. U.
,
2017
, “
Experimental Characterization of Stable Liquid Rivulets on Inclined Surfaces: Influence of Surface Tension, Viscosity and Inclination Angle on the Interfacial Area
,”
Chem. Eng. Res. Des.
,
125
, pp.
226
232
.10.1016/j.cherd.2017.07.022
16.
Singh
,
R. K.
,
Galvin
,
J. E.
, and
Sun
,
X.
,
2016
, “
Three-Dimensional Simulation of Rivulet and Film Flows Over an Inclined Plate: Effects of Solvent Properties and Contact Angle
,”
Chem. Eng. Sci.
,
142
, pp.
244
257
.10.1016/j.ces.2015.11.029
17.
Singh
,
R. K.
,
Galvin
,
J.
,
Whyatt
,
G. A.
, and
Sun
,
X.
,
2017
, “
Breakup of a Liquid Rivulet Falling Over an Inclined Plate: Identification of a Critical Weber Number
,”
Phys. Fluids.
, 29(5), p. 052101.10.1063/1.4981920
18.
Meredith
,
K. V.
,
Heather
,
A.
,
de Vries
,
J.
, and
Xin
,
Y.
,
2011
, “
A Numerical Model for Partially-Wetted Flow of Thin Liquid Film
,”
WIT Trans. Eng. Sci.
,
70
, pp.
239
250
.10.2495/MPF110201
19.
OpenFoam,
2018
,
The Open Source CFD Toolbox, (n.d.)
,”
OpenFoam, London, UK, accessed Mar. 10, 2019,
http://www.openfoam.org
20.
Martin
,
M.
,
Defraeye
,
T.
,
Derome
,
D.
, and
Carmeliet
,
J.
,
2015
, “
A Film Flow Model for Analysing Gravity-Driven, Thin Wavy Fluid Films
,”
Int. J. Multiph. Flow
,
73
, pp.
207
216
.10.1016/j.ijmultiphaseflow.2015.03.010
21.
Singh
,
K.
,
Sharabi
,
M.
,
Ambrose
,
S.
,
Eastwick
,
C.
, and
Jefferson-Loveday
,
R.
,
2019
, “
Prediction of Film Thickness of an Aero-Engine Bearing Chamber Using Coupled VOF and Thin Film Model
,”
ASME
Paper No. GT2019-91314. 10.1115/GT2019-91314
22.
Singh
,
K.
,
Sharabi
,
M.
,
Ambrose
,
S.
,
Eastwick
,
C.
,
Jefferson-Loveday
,
R.
,
Cao
,
J.
, and
Jacobs
,
A.
,
2019
, “
Assessment of an Enhanced Thin Film Model to Capture Wetting and Drying Behavior in an Aero-Engine Bearing Chamber
,”
ASME
Paper No. GT2019-91323. 10.1115/GT2019-91323
23.
Bristot
,
A.
,
2017
,
Application of the Volume of Fluid Method With Heat Transfer to a Two- Shaft Aero-Engine Bearing Chamber
,
The University of Nottingham
,
Nottingham, UK
24.
Kakimpa
,
B.
,
Morvan
,
H.
, and
Hibberd
,
P. S.
,
2015
, “
Solution Strategies for Thin Film Rimming Flow Modelling
,”
ASME
Paper No. GT2015-43503. 10.1115/GT2015-43503
25.
ANSYS (US)
,
2020
, “
ANSYS Fluent Theory Guide, release2020R1
,” ANSYS, Canonsburg, PA.
26.
Chandra
,
B.
,
Collicott
,
S. H.
, and
Munson
,
J. H.
,
2013
, “
Scavenge Flow in a Bearing Chamber With Tangential Sump Off-Take
,”
ASME J. Eng. Gas Turbines Power
,
135
(
3
), p.
032503
.10.1115/1.4007869
27.
Kakimpa
,
B.
,
Morvan
,
H.
, and
Hibberd
,
S.
,
2016
, “
The Depth-Averaged Numerical Simulation of Laminar Thin-Film Flows With Capillary Waves
,”
ASME J. Eng. Gas Turbines Power
,
138
(
11
), p.
112501
.10.1115/1.4033471
28.
Moyle
,
D. T.
,
Chen
,
M. S.
, and
Homsy
,
G. M.
,
1999
, “
Nonlinear Rivulet Dynamics During Unstable Wetting Flows
,”
Int. J. Multiphas Flow
,
25
(
6–7
), pp.
1243
1262
.10.1016/S0301-9322(99)00062-2
29.
Lukyanov
,
A. V.
, and
Shikhmurzaev
,
Y. D.
,
2007
, “
Effect of Flow Field and Geometry on the Dynamic Contact Angle
,”
Phys. Rev. E. Stat. Nonlin. Soft Matter Phys
,
75
(
5
), p.
51604
.10.1103/PhysRevE.75.051604
30.
Takamasa
,
T.
, and
Hazuku
,
T.
,
2000
, “
Measuring Interfacial Waves on Film Flowing Down a Vertical Plate Wall in the Entry Region Using Laser Focus Displacement Meters
,”
Int. J. Heat Mass Transfer
,
43
(
15
), pp.
2807
2819
.10.1016/S0017-9310(99)00335-X
31.
Patankar
,
S. V.
,
1990
,
Numerical Heat Transfer and Fluid Flow
,
Hemisphere Publishing Corporation
,
Washington, DC
.
32.
Celik
,
I. B.
, and
Li
,
J.
,
2005
, “
Assessment of Numerical Uncertainty for the Calculations of Turbulent Flow Over a Backward-Facing Step
,”
Int. J. Numer. Methods Fluids
,
49
(
9
), pp.
1015
1031
.10.1002/fld.1040
33.
Kurz
,
W.
,
Dullenkopf
,
K.
, and
Bauer
,
H.-J.
,
2013
, “
Capacitive Film Thickness Measurements in a Ventless Aero-Engine Bearing Chamber—Influence of Operating Conditions and Offtake Design
,”
ASME J. Eng. Gas Turbines Power
,
135
(
11
), p.
112504
.10.1115/1.4025067