Modeling of Partially Wetting Liquid Film Using an Enhanced Thin Film Model for Aero-Engine Bearing Chamber Applications

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 [1–4]. 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 1 mm, 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.

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.

where ρ_l is the liquid density, h is the film thickness, $∇s$ is the surface gradient operator, u_l is the velocity of the liquid, and S_m 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.

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, S_i,mom, are employed in this study.

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.

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.

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.

where ρ_l is the liquid density, μ_l is the liquid viscosity, ν_l is the liquid kinematic viscosity, g_x is the gravitational acceleration, and q_l 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].

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.

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⁻⁵ for the momentum and continuity equations.

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.

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 (p_avg.) 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 (p_avg.). 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.7–2.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.

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%.

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.

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].

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.

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.

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.

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.

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.

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/m³, σ = 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 ß.

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.

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.

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.

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.

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.

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.

• 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).

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, m³/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

β =

empirical parameter to account for contact angle force

δ =

standard deviation

ε21 =

h2–h1

θ =

contact angle, deg

μ =

dynamic viscosity, N-s/m²

ν =

kinematic viscosity, m²/s

ρ =

density, kg/m³

σ =

surface tension, N/m

τ =

shear stress, N/m²

ϕ =

inclination angle, deg

Subscripts

c =

critical

f =

film

l =

liquid

m =

mean

s =

starting

Acronyms

CD =

computational domain

ETFM =

Eulerian thin/wall film model

LPM =

liter per minute

VOF =

volume of fluid