*3.1. Governing Equations*

In the present study, the flow field transits from single-phase flow to air–oil two-phase flow with the increasing speed of the grooved disk. Thus, a three-dimensional CFD approach dealing with a multiphase flow problem presented for the open grooved two-disk system [31] was utilized in this flow field simulation. In the visualization experiment, it was found that the air–oil interface between the inner and the outer radii of the grooved disk was rather stable after the rotating speed of the grooved disk was reached and kept constant. Thus, a steady-state model was established. The model was developed using the commercial software ANSYS FLUENT 15.0. In this model, the VOF approach developed by Hirt and Nichols [38] was adopted to track the two-phase interface of the air and the oil. A volume fraction ϕoil was used to mark the volume fraction of the oil phase. In the model, ϕoil = 1 represents the pure oil phase in a cell, while ϕoil = 0 represents the cell full of air phase. If 0 < ϕoil < 1, it describes the air–oil two-phase interface. The volume fractions of two phases satisfy the following constraint:

$$
\varphi\_{\dot{a}\dot{r}} + \varphi\_{o\ddot{o}l} = 1 \tag{1}
$$

In the VOF method, the flow properties are averaged by phase volume fraction [39]. The density, the dynamic viscosity, and the effective thermal conductivity coefficients of the air–oil two-phase flow are given by:

$$
\rho = \varphi\_{\rm oil} \rho\_{\rm oil} + (1 - \varphi\_{\rm oil}) \rho\_{\rm air} \tag{2}
$$

$$
\mu = \varphi\_{a\text{il}}\mu\_{oil} + (1 - \varphi\_{a\text{il}})\mu\_{air} \tag{3}
$$

$$k\_{eff} = \varphi\_{oil} k\_{oil} + (1 - \varphi\_{oil}) k\_{air} \tag{4}$$

The governing equations of the CFD method are the momentum equation, the continuity equation, and the energy equation. These equations are solved in the VOF multiphase model. The VOF model solves one single momentum equation, and the velocity field results are shared among both phases and result in a direct velocity coupling of the phases at the interface [34]. This computational model selects steady reference frame. The simulation is steady-state, and the transient terms in the governing equations are omitted.

The continuity equation is:

$$\nabla \cdot \left( \rho \overrightarrow{\boldsymbol{\nu}} \right) = 0 \tag{5}$$

The momentum equation is:

$$\nabla \cdot \left( \rho \stackrel{\scriptstyle \rightarrow}{\boldsymbol{\nu}} \stackrel{\scriptstyle \rightarrow}{\boldsymbol{\nu}} \right) = \nabla \cdot \left( \mu \left( \boldsymbol{\nabla} \stackrel{\scriptstyle \rightarrow}{\boldsymbol{\nu}} + \left( \boldsymbol{\nabla} \stackrel{\scriptstyle \rightarrow}{\boldsymbol{\nu}} \right)^{T} \right) \right) - \nabla \cdot \boldsymbol{p} + \rho \stackrel{\scriptstyle \rightarrow}{\boldsymbol{g}} + \stackrel{\scriptstyle \rightarrow}{F} \tag{6}$$

Equations (2) and (3) were implemented as the averaged dynamic viscosity and the density of the fluid in the above equations. <sup>→</sup> υ denotes the velocity vector of fluid, *p* is the pressure of fluid, <sup>→</sup> *g* is the gravity acceleration vector, and <sup>→</sup> *F* denotes the surface tension source term. The surface tension force is given as:

$$\overrightarrow{F\_i} = \sigma \frac{\rho \mathbf{x}\_i \nabla \varphi\_i}{(\rho\_i + \rho\_j)/2} \tag{7}$$

where *i* and *j* denote the air and the oil phase, respectively. σ is the surface tension coefficient assumed to be constant 0.03 N/m in this simulation. κ*<sup>i</sup>* is the surface curvature at the interface defined as the divergence of the unit normal. For two phases presented in a cell, κ*<sup>i</sup>* = −κ*j*.

$$\kappa\_i = \nabla \cdot \hbar\_i \tag{8}$$

where:

$$
\hbar\_i = \hbar\_w \cos \theta\_w + \hat{t}\_w \sin \theta\_w \tag{9}
$$

θ*<sup>w</sup>* is contact angle at the wall. θ*<sup>w</sup>* = 9.3◦ and 10.1◦ on the surfaces of the stationary disk made of quartz glass and the rotating aluminum disk, respectively. The data were obtained from the contact angle test results described in reference [40]. *n*ˆ *<sup>w</sup>* and ˆ*tw* are the unit vectors normal and tangential to the interface, respectively. Since the surface tension coefficient was assumed to be constant, the Marangoni effect was neglected due to zero surface tension gradient.

In the simulation, the air phase was set as the primary phase. The continuity equation of the air phase was solved first, given as:

$$\nabla \cdot \left( \alpha\_{\dot{a}\dot{r}} \rho\_{\dot{a}\dot{r}} \stackrel{\rightarrow}{\upsilon}\_{a\dot{r}} \right) = 0 \tag{10}$$

After obtaining the volume fraction of the air phase in a cell, the volume fraction of the oil was determined from the constraint (1) as:

$$
\rho\_{oil} = 1 - \rho\_{air} \tag{11}
$$

The energy equation omitting radiation as a source term for the fluid domain was deduced as [41]:

$$\nabla \cdot (\overrightarrow{\boldsymbol{\nu}} \, \boldsymbol{\rho} \, \boldsymbol{E}) = \nabla \cdot (\boldsymbol{k}\_{\boldsymbol{\varepsilon} \, \boldsymbol{f} \, \boldsymbol{\Gamma}} \, \boldsymbol{\Gamma} \, \boldsymbol{T}) \tag{12}$$

where:

$$E = \frac{\rho\_{\rm oil}\rho\_{\rm oil}E\_{\rm oil} + (1 - \rho\_{\rm oil})\rho\_{\rm air}E\_{\rm air}}{\rho\_{\rm oil}\rho\_{\rm oil} + (1 - \rho\_{\rm oil})\rho\_{\rm air}} \tag{13}$$

The energy equation for the solid domain is:

$$\nabla \cdot (k\_s \nabla T) = 0 \tag{14}$$

where *T* denotes the temperature shared between the air and the oil phases, *ke*ff represents the effective thermal conductivity of the two-phase flow given by Equation (4), ϕ is the fluid volume fraction, *E* is the specific sensible enthalpy, and *ks* denotes the thermal conductivity of the solid structure. The oil and air subscripts identify the oil and the air parameters, respectively.
