Analytical and numerical validation of a plate–plate tribometer for measuring wall slip

Abstract

We model the slip length tribometer (SLT), originally presented by Pelz et al. (J. Fluid Mech., vol. 948, 2022, p. A8) in OpenFOAM. The plate tribometer is especially designed to simultaneously measure viscosity and slip length for lubrication gaps in the range of approximately 10 $\mathrm {\mu }$m at temperatures and surface roughnesses relevant to technical applications, with a temperature range of $-30$ to $100\,^\circ \mathrm {C}$ and surface roughness ranging from $10\ \mathrm {nm}$ to $1\ \mathrm {\mu }\mathrm {m}$. A simplified analytical model presented by Pelz et al. (J. Fluid Mech., vol. 948, 2022, p. A8) infers the slip length of the plate from the experimentally measured torque and the plate gap height. The present work verifies the analytical model using axisymmetric flow simulations and presents the effect of inlet on the numerical velocity profiles. The simulation results are in very good agreement with the results of the analytical model. The main conclusion drawn from this study is the validation of the Navier-slip boundary condition as an effective model for partial slip in computational fluid dynamics simulations and the negligible influence of the inlet on the fluid flow between the SLT's plates.

Impact Statement

This research explores the effect of partial slip on the flow of fluids in a microfluidic device called the slip length tribometer (SLT). By using a simplified version of Navier–Stokes equations with the Navier-slip boundary condition, partial slip on the SLT's plates was successfully correlated with experimental measurements. We use the OpenFOAM open-source computational fluid dynamics software to perform more detailed simulations that include the effect of a flow inlet and an outflow at the outer radius between the SLT's plates. The more detailed simulation confirms the experimental measurement and the simplified model and provides a basis for future research on partial slip caused by surface roughness.

1. Introduction

For many flow problems in science or engineering, the no-slip boundary condition is applied to model the interaction between fluid flow and surrounding solid walls (e.g. Spurk & Aksel 2019). The no-slip condition states that the fluid velocity at an impermeable wall is identical to the velocity of the wall itself, preventing relative motion between the wall and the fluid molecules at the wall. Even if the no-slip boundary condition is an appropriate model for many technical applications, deviations have been observed for microfluidic flows (cf. Lauga, Brenner & Stone 2007; Neto et al. 2005) or dynamic wetting flows, where the breakdown of the no-slip boundary condition has been clearly pointed out by Huh & Scriven (1971). This insight is by no means new, as the slip boundary condition was already postulated by Navier (1822). He postulated that the fluid slips, with a relative velocity being linearly proportional to the shear rate at the wall. Here, the constant of proportionality is the slip length, which, according to Helmholtz & von Piotrowski (1860), can be interpreted as an effective increase in gap clearance. The Navier-slip boundary condition expresses a balance between the friction force (opposing the tangential motion) and the component of the viscous shear force parallel to the solid wall, i.e.

(1.1)\begin{equation} -\beta ( \boldsymbol{v} - \boldsymbol{v}_w )_\parallel{ = } \boldsymbol{(Sn)}_{||}. \end{equation}

Here, $\beta$ is a coefficient describing the amount of friction between liquid and solid particles, $\boldsymbol {S} = \mu (\boldsymbol {\nabla } \boldsymbol {v} + \boldsymbol {\nabla } \boldsymbol {v}^{\textsf {T}})$ is the viscous stress tensor and $\mu$ denotes the dynamic viscosity. Within continuum thermodynamics, (1.1) is the simplest, namely linear, closure to model a priori the unknown relative tangential velocity which appears in the entropy production contribution due to relative (tangential) motion between two phases; see Bothe (2022) for more details. Moreover, the solid wall is assumed to be impermeable, which implies that the normal component of the relative velocity vanishes, i.e.

(1.2)\begin{equation} ( \boldsymbol{v} - \boldsymbol{v}_w ) \boldsymbol{\cdot} \boldsymbol{n} = 0, \end{equation}

where $\boldsymbol {n}$ denotes the outer normal vector field to the solid wall. Note that the full boundary condition for solving the Navier–Stokes (or Stokes) equations for fluid flow in a confined domain with Navier slip is given by (1.1) and (1.2) together.

An equivalent formulation of (1.1) is obtained by dividing the equation by the friction coefficient, leading to

(1.3)\begin{equation} ( \boldsymbol{v} - \boldsymbol{v}_w )_\parallel{ + } 2 \lambda \boldsymbol{(Dn)}_{||} = 0. \end{equation}

Here, $\boldsymbol {D}=\frac 12 (\boldsymbol {\nabla } \boldsymbol {v}+ (\boldsymbol {\nabla } \boldsymbol {v})^{\textsf {T}})$ is the rate-of-deformation tensor and the quantity

(1.4)\begin{equation} \lambda : = \frac{\mu}{\beta} \end{equation}

is called the slip length. It is known from molecular dynamics simulations (cf. Mehrnia & Pelz 2021) and experimental investigation that the value of the slip length is typically of the order of 10 to 100 nm (see Zhu & Granick 2002; Bonaccurso, Butt & Craig 2003; Granick, Zhu & Lee 2003; Neto et al. 2005; Vinogradova & Yakubov 2006). Hence, the slip length is far below the characteristic length scales of the flow for many technical applications. Furthermore, the impact of $\lambda$ on the macroscopic flow may be small, particularly for single-phase flows. However, for dynamic wetting flows or flow processes where the characteristic length scale is of the order of the slip length, e.g. sealing systems and hydraulic components, it is important to determine the value of the slip parameter. More details of the latest research can be found in the review articles by Neto et al. (2005), Cao et al. (2009), Vinogradova & Belyaev (2011) and Pelz et al. (2022). Here, the slip parameter is of the order of magnitude of the prescribing surface roughness. It is stated, especially by Pelz et al. (2022), that most of the experimental techniques used to measure the slip parameter are inadequate for determining slip length for technically relevant rough surfaces. This is mainly due to the small characteristic length scale of the measurement device, being of the order of 10 $\mathrm {\mu }$m. To adequately measure the slip length for large surface roughnesses, a much higher characteristic length scale of the system is needed.

For this reason, the slip length tribometer (SLT) was developed (see Pelz et al. (2022) and Corneli (2022) for details). The SLT is a plate–plate tribometer designed to measure both the lubrication gap height and the frictional torque transmitted from the rotating plate, through the entrapped liquid, to the stationary plate. It overcomes limitations of existing devices like the surface force apparatus (Horn et al. 2000) and the atomic force microscope (Henry & Craig 2009), particularly in characterizing technically relevant rough surfaces at varying temperatures. The SLT can analyse machined surfaces with arithmetic mean roughnesses ranging from $10 \ \mathrm {nm}$ to $1 \ \mathrm {\mu }\mathrm {m}$, across temperatures from $-30$ to $100\,^\circ \mathrm {C}$. Its ability to handle such large roughnesses is largely due to the significant characteristic length scale, i.e. the diameter of the plates, being of the order of 10 mm. This results in a wetted area of approximately $3 \times 10^3\ \mathrm {mm}^2$, which is six orders of magnitude greater than that of existing apparatus for wall slip measurements, enabling an in situ reduction of uncertainty. Additional reduction of overall uncertainty is achieved through in situ calibration of the capacitive position sensors using an interferometric white-light thickness measurement device. Furthermore, stochastic uncertainty is minimized through an automated measurement process, utilizing between 300 and 800 measurement points, where each point is averaged over 70 000 individual measurements within a 10 s interval for a single slip length measurement.

The slip length is derived by using an analytical solution of Navier–Stokes equations with Navier-slip boundary condition for a Couette-type flow between two rotating plates, cf. (1.5), relating the measured torque and gap height to the slip length itself.

However, the simplified approach, derived from a solution of the Stokes equation, neglects the feed flow entering and leaving the device during measurement. Therefore, the main goal of the present work is to show, utilizing computational fluid dynamics simulations in the open-source software OpenFOAM-v2312 (2023) (see Marić, Höpken & Mooney (2014) for workings and usage of OpenFOAM), that the simplified analytical model

(1.5)\begin{equation} M = \frac{\mu I_p \varOmega}{h + 2\lambda} \end{equation}

for the torque $M$ as a function of the plate distance $h$, the polar moment of the area $I_p = {{\rm \pi} }R^4/2$ of the lower plate and angular velocity $\varOmega$ is sufficient to infer the slip length $\lambda$ from the experimentally measured torque. The structure of the paper is as follows. In § 2, the simplified analytical model for the torque is presented, followed by a description of the simulation models and set-up. Section 3 discusses the results and provides a comparative analysis of the employed models. The presented work is then summarized in § 4.

2. Materials and methods

2.1 Analytical model

We consider the steady and rotationally symmetric flow of an incompressible Newtonian fluid between two plates of radius R. The upper plate rotates with constant angular velocity $\varOmega$. The lower plate is at rest. Let $(r, \phi, z)$ be the cylindrical coordinates so that the plates are at $z = 0$ and $z = h$. Figure 1 shows a sketch of the tribometer.

Figure 1. Sketch of the principle of the analytical model (see Corneli (2022) for details).

2.1.1 Equations of motion

Compared with the work of Corneli (2022), where a simpler linear model is assumed, we apply Navier-slip boundary conditions to full steady-state Navier–Stokes equations in cylindrical coordinates.

We introduce the velocity vector $\boldsymbol {v}=v_r \boldsymbol {e_r}+v_{\phi } \boldsymbol {e_{\phi }}+v_z \boldsymbol {e_z}$ and the pressure $p$. Considering the rotational symmetry, the Navier–Stokes equations have the following form (Spurk & Aksel 2019):

(2.1)\begin{gather} \frac{\partial v_r}{\partial r} + \frac{\partial v_z}{\partial z} + \frac{v_r}{r} = 0, \end{gather}

(2.2)\begin{gather}v_r \frac{\partial v_r}{\partial r} + v_z \frac{\partial v_r}{\partial z} - \frac{v_{\phi}^2}{r} = {-}\frac{1}{\rho} \frac{\partial p}{\partial r} + \nu \left(\frac{\partial^2 v_r}{\partial r^2} + \frac{1}{r}\frac{\partial v_r}{\partial r} - \frac{v_r}{r^2} + \frac{\partial^2 v_r}{\partial z^2} \right), \end{gather}

(2.3)\begin{gather}v_r \frac{\partial v_{\phi}}{\partial r} + v_z \frac{\partial v_{\phi}}{\partial z} + \frac{v_r v_{\phi}}{r} = \nu \left(\frac{\partial^2 v_{\phi}}{\partial r^2} + \frac{1}{r}\frac{\partial v_{\phi}}{\partial r} - \frac{v_{\phi}}{r^2} + \frac{\partial^2 v_{\phi}}{\partial z^2} \right), \end{gather}

(2.4)\begin{gather}v_r \frac{\partial v_z}{\partial r} + v_z \frac{\partial v_z}{\partial z} = {-}\frac{1}{\rho} \frac{\partial p}{\partial z} + \nu \left(\frac{\partial^2 v_z}{\partial r^2} + \frac{1}{r}\frac{\partial v_z}{\partial r} + \frac{\partial^2 v_z}{\partial z^2} \right). \end{gather}

Here, $\rho$ denotes the density and $\nu$ is the kinematic viscosity.

2.1.2 Non-dimensionalization of the equations of motion

The plate radius $R$ and the gap width $h$ are the characteristic lengths used to non-dimensionalize the equations. We define the dimensionless variables

(2.5a–f)\begin{equation} \hat{r} : = \frac{r}{R}, \quad \hat{z} : = \frac{z}{h}, \quad \hat{v}_r : = \frac{v_r}{R \varOmega}, \quad \hat{v}_{\phi} : = \frac{v_{\phi}}{R \varOmega}, \quad \hat{v}_z: = \frac{v_z}{h \varOmega}, \quad \hat{p}: = \frac{p}{\rho R^2 \varOmega^2}. \end{equation}

Substituting the dimensionless variables into the dimensioned equations (2.1)–(2.4) yields the dimensionless system of equations as

(2.6)\begin{gather} \frac{\partial \hat{v}_r}{\partial \hat{r}} + \frac{\partial \hat{v}_z}{\partial \hat{z}} + \frac{\hat{v}_r}{\hat{r}} = 0, \end{gather}

(2.7) \begin{gather}\hat{v}_r \frac{\partial \hat{v}_r}{\partial \hat{r}} + \hat{v}_z \frac{\partial \hat{v}_z}{\partial \hat{z}} - \frac{\hat{v}_{\phi}^2}{\hat{r}} = {-} \frac{\partial \hat{p}}{\partial \hat{r}} + \frac{1}{Re} \left[ \delta^2 \left( \frac{\partial^2 \hat{v}_r}{\partial \hat{r}^2} + \frac{1}{\hat{r}}\frac{\partial \hat{v}_r}{\partial \hat{r}} - \frac{\hat{v}_r}{\hat{r}^2} \right) + \frac{\partial^2 \hat{v}_r}{\partial \hat{z}^2} \right], \end{gather}

(2.8)\begin{gather}\hat{v}_r \frac{\partial \hat{v}_{\phi}}{\partial \hat{r}} + \hat{v}_z \frac{\partial \hat{v}_{\phi}}{\partial \hat{z}} + \frac{\hat{v}_r \hat{v}_{\phi}}{\hat{r}} = \frac{1}{Re} \left[ \delta^2 \left( \frac{\partial^2 \hat{v}_{\phi}}{\partial \hat{r}^2} + \frac{1}{\hat{r}}\frac{\partial \hat{v}_{\phi}}{\partial \hat{r}} - \frac{\hat{v}_{\phi}}{\hat{r}^2} \right) + \frac{\partial^2 \hat{v}_{\phi}}{\partial \hat{z}^2} \right], \end{gather}

(2.9)\begin{gather}\hat{v}_r \frac{\partial \hat{v}_z}{\partial \hat{r}} + \hat{v}_z \frac{\partial \hat{v}_z}{\partial \hat{z}} = {-} \frac{1}{\delta^2} \frac{\partial \hat{p}}{\partial \hat{z}} + \frac{1}{Re} \left[ \delta^2 \left( \frac{\partial^2 \hat{v}_z}{\partial \hat{r}^2} + \frac{1}{r}\frac{\partial \hat{v}_z}{\partial \hat{r}}\right) + \frac{\partial^2 \hat{v}_z}{\partial \hat{z}^2} \right]. \end{gather}

The non-dimensional parameters used in this study are the Reynolds number $Re=({h^2 \varOmega })/{\nu }$ and the plate spacing $\delta ={h}/{R}$. According to the operation parameters of the SLT, the Reynolds number is of the order of $10^{-9}$. Therefore, the inertial terms can be neglected, and the following Stokes equations apply:

(2.10)\begin{gather} \frac{\partial \hat{v}_r}{\partial \hat{r}} + \frac{\partial \hat{v}_z}{\partial \hat{z}} + \frac{\hat{v}_r}{\hat{r}} = 0, \end{gather}

(2.11)\begin{gather}- \frac{\partial \hat{p}}{\partial \hat{r}} + \frac{1}{Re} \left[ \delta^2 \left( \frac{\partial^2 \hat{v}_r}{\partial \hat{r}^2} + \frac{1}{\hat{r}}\frac{\partial \hat{v}_r}{\partial \hat{r}} - \frac{\hat{v}_r}{\hat{r}^2} \right) + \frac{\partial^2 \hat{v}_r}{\partial \hat{z}^2} \right] = 0, \end{gather}

(2.12)\begin{gather}\frac{1}{Re} \left[ \delta^2 \left( \frac{\partial^2 \hat{v}_{\phi}}{\partial \hat{r}^2} + \frac{1}{\hat{r}}\frac{\partial \hat{v}_{\phi}}{\partial \hat{r}} - \frac{\hat{v}_{\phi}}{\hat{r}^2} \right) + \frac{\partial^2 \hat{v}_{\phi}}{\partial \hat{z}^2} \right] = 0, \end{gather}

(2.13)\begin{gather}- \frac{1}{\delta^2} \frac{\partial \hat{p}}{\partial \hat{z}} + \frac{1}{Re} \left[ \delta^2 \left( \frac{\partial^2 \hat{v}_z}{\partial \hat{r}^2} + \frac{1}{\hat{r}}\frac{\partial \hat{v}_z}{\partial \hat{r}}\right) + \frac{\partial^2 \hat{v}_z}{\partial \hat{z}^2} \right] = 0. \end{gather}

It should be emphasized that (2.12) is uncoupled and can be solved independently. Since we are only interested in the torque at the bottom plate resulting from the shear stress $\tau _{\phi z}=\mu ({\partial v_{\phi }}/{\partial z})$, only the circumferential velocity component $\hat {v}_\phi (\hat {r},\hat {z})$ will be computed. For simplicity, the circumferential velocity component $\hat {v}_\phi (\hat {r},\hat {z})$ is denoted as $\hat {v}(\hat {r},\hat {z})$. The equation

(2.14)\begin{equation} \delta^2 \left( \frac{\partial^2 \hat{v}}{\partial \hat{r}^2} + \frac{1}{\hat{r}}\frac{\partial \hat{v}}{\partial \hat{r}} - \frac{\hat{v}}{\hat{r}^2} \right) + \frac{\partial^2 \hat{v}}{\partial \hat{z}^2} = 0\end{equation}

is solved in the following.

2.1.3 Boundary conditions

In order to solve the equation of motion (2.14) uniquely, the boundary conditions must be specified. The Navier-slip boundary condition is applied at the plates, and considering the rotational symmetry, the following boundary conditions apply:

(2.15)\begin{gather} \hat{v}(\hat{r} = 0,\hat{z}) = 0, \end{gather}

(2.16)\begin{gather}- \left.\frac{\lambda_u}{h} \frac{\partial \hat{v}(\hat{r},\hat{z})}{\partial \hat{z}} \right|_{\hat{z} = 0} + \hat{v}(\hat{r},\hat{z} = 0) = 0 \end{gather}

and

(2.17)\begin{equation} \left.\frac{\lambda_o}{h} \frac{\partial \hat{v}(\hat{r},\hat{z})}{\partial \hat{z}} \right|_{\hat{z} = 1} + \hat{v}(\hat{r},\hat{z} = 1) = \hat{r}. \end{equation}

2.1.4 Analytical solution

Substituting $\hat {v}=\hat {r} f(\hat {z})$ in the differential equation (2.14), with which the symmetry condition (2.15) is automatically satisfied, leads to

(2.18)\begin{equation} \frac{\partial^2 f(\hat{z})}{\partial \hat{z}^2} = 0.\end{equation}

The boundary conditions (2.16) and (2.17) simplify to

(2.19)\begin{equation} - \left.\frac{\lambda_u}{h} \frac{\partial \hat{f}(\hat{z})}{\partial \hat{z}} \right|_{\hat{z} = 0} + \hat{f}(0) = 0 \end{equation}

and

(2.20)\begin{equation} \left.\frac{\lambda_o}{h} \frac{\partial \hat{f}(\hat{z})}{\partial \hat{z}} \right|_{\hat{z} = 1} + \hat{f}(1) = 1. \end{equation}

The integration of the differential equation (2.18) considering the boundary conditions (2.19) and (2.20) leads to the velocity profile

(2.21)\begin{equation} \hat{v}(\hat{r},\hat{z}) = \frac{\hat{r}}{1 + (\lambda_o + \lambda_u)/h}\left(\hat{z} + \frac{\lambda_u}{h} \right). \end{equation}

Since the fluid is Newtonian, the shear stress $\tau _{\phi z}$ (see Spurk & Aksel (2019) for details) can be expressed as

(2.22)\begin{equation} \tau_{\phi z} = \mu \frac{\partial v_{\phi}}{\partial z} = \frac{\mu \varOmega r}{h + \lambda_u + \lambda_o}. \end{equation}

By integrating the shear stress $\tau _{\phi z}$ over the area of the bottom plate, the torque becomes

(2.23)\begin{equation} M = \int_{\phi = 0}^{\phi = 2{\rm \pi}} \int_{r = 0}^{r = R} r\tau_{\phi z} r \, \mathrm{d}r \, \mathrm{d}\phi = \frac{\mu I_p \varOmega}{h + \lambda_u + \lambda_o}.\end{equation}

For this study, we consider a homogeneous material for both rheometer plates; hence, the slip length is assumed to be equal for both plates. In this case, the formula for the torque simplifies to (1.5).

For a given measured torque at a specific gap height $h$, the linear relationship (1.5) between $M^{-1}$ and $h$ provides the value of slip length (see Kuhr, Corneli & Pelz (2022) and Pelz et al. (2022) for details).

2.2 Simulation model

Figure 2(a) shows a sketch of the principle of the SLT, simulated using the open-source software (OpenFOAM-v2312 2023). Two simulation models are presented: the actual model and the simplified model. The actual model has a cylindrical tube that drives the fluid between the plates, through which the fluid (density $\rho$, dynamic viscosity $\mu$) flows in with constant volumetric flow rate $Q$ at pressure $p_1$. The actual simulation model is used to estimate the effect of the cylindrical tube on the resulting torque. The cylindrical tube is modelled as a hole (radius $R_1$) in the simplified simulation model and the inlet is present at the top rheometer plate. A schematic diagram of the axisymmetric computational domain is shown in figure 2(b). At the outlet, the fluid leaves the computational domain. At the solid boundary, the simulation utilizes the OpenFOAM implementation of Navier slip developed by Gründing, Raju & Maric (2022) with the slip length provided by the analytical model.

Figure 2. (a) Sketch of the principle of the simulation model (see Corneli (2022) for details). (b) Schematic diagram of an axisymmetric domain of the actual model with the tribometer radius $R$, gap height $h$, an inlet height $h_1$ and inlet radius $R_1$. For the simplified simulation model, the inlet is present at the top tribometer plate, i.e. the blue region is excluded from the computational domain.

2.2.1 Equations of motion

Since OpenFOAM does not support the formulation of the balance equations in cylindrical coordinates, we consider the Navier–Stokes equations in Cartesian coordinates.

For the velocity vector $\boldsymbol {v}=[v_x, v_y, v_z]^{\textsf {T}})$ and the pressure $p$, the following balance equations for mass and momentum (Spurk & Aksel 2019) apply in the computational domain:

(2.24)\begin{gather} \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} + \frac{\partial v_z}{\partial z} = 0, \end{gather}

(2.25) \begin{gather} \left.\begin{array}{c} \displaystyle \rho \left( v_x \dfrac{\partial v_x}{\partial x} + v_y \dfrac{\partial v_x}{\partial y} + v_z \dfrac{\partial v_x}{\partial z} \right) = {-}\dfrac{\partial p}{\partial x} + \mu \left( \dfrac{\partial^2 v_x}{\partial x^2} + \dfrac{\partial^2 v_x}{\partial y^2} + \dfrac{\partial^2 v_x}{\partial z^2} \right), \\ \displaystyle \rho \left( v_x \dfrac{\partial v_y}{\partial x} + v_y \dfrac{\partial v_y}{\partial y} + v_z \dfrac{\partial v_y}{\partial z} \right) = {-}\dfrac{\partial p}{\partial y} + \mu \left( \dfrac{\partial^2 v_y}{\partial x^2} + \dfrac{\partial^2 v_y}{\partial y^2} + \dfrac{\partial^2 v_y}{\partial z^2} \right), \\ \displaystyle \rho \left( v_x \dfrac{\partial v_z}{\partial x} + v_y \dfrac{\partial v_z}{\partial y} + v_z \dfrac{\partial v_z}{\partial z} \right) = {-}\dfrac{\partial p}{\partial z} + \mu \left( \dfrac{\partial^2 v_z}{\partial x^2} + \dfrac{\partial^2 v_z}{\partial y^2} + \dfrac{\partial^2 v_z}{\partial z^2} \right). \end{array}\right\} \end{gather}

2.2.2 Boundary conditions

The Navier-slip boundary condition and the condition of impermeability hold at the bottom and top plates and the inlet channel wall, and they can then be expressed in the form

(2.26) \begin{gather} - \lambda_u \left(\begin{array}{c} \displaystyle \dfrac{\partial v_x}{\partial z} + \dfrac{\partial v_z}{\partial x} \\ \displaystyle \dfrac{\partial v_y}{\partial z} + \dfrac{\partial v_z}{\partial y} \\ 0 \end{array}\right)_{z = 0} + \left(\begin{array}{c} v_x \\ v_y \\ v_z \end{array}\right)_{z = 0} = \left(\begin{array}{c} 0 \\ 0 \\ 0\end{array}\right), \end{gather}

(2.27) \begin{gather}\lambda_o \left(\begin{array}{c} \displaystyle \dfrac{\partial v_x}{\partial z} + \dfrac{\partial v_z}{\partial x} \\ \displaystyle \dfrac{\partial v_y}{\partial z} + \dfrac{\partial v_z}{\partial y} \\ 0 \end{array}\right)_{z = h} + \left(\begin{array}{c} v_x \\ v_y \\ v_z \end{array}\right)_{z = h} = \left(\begin{array}{c} 0 \\ r \varOmega \\ 0\end{array}\right)_{z = h}, \end{gather}

(2.28) \begin{gather}- \lambda_i \left(\begin{array}{c} 0 \\ \displaystyle \dfrac{\partial v_x}{\partial y} + \dfrac{\partial v_y}{\partial x} \\ \displaystyle \dfrac{\partial v_x}{\partial z} + \dfrac{\partial v_z}{\partial x} \end{array}\right)_{x = R_1} + \left(\begin{array}{c} v_x \\ v_y \\ v_z \end{array}\right)_{x = R_1} = \left(\begin{array}{c} 0 \\ 0 \\ 0\end{array}\right)_{x = R_1}, \end{gather}

where $r=\sqrt {x^2+y^2}$ describes the distance between any point $P$ on the upper plate with coordinates $(x,y,h)$ and the axis of rotation of the SLT. At the inlet, a constant pressure of value $p_i$ is applied as a boundary condition, i.e.

(2.29)\begin{equation} p(x,y,h + h_1) = p_i \quad \text{for} \ 0 \leq \sqrt{x^2 + y^2} \leq R_1, \end{equation}

where Corneli (2022) relates $p_i$ to the plate height $h$ as

(2.30)\begin{equation} h = \frac{\rm \pi}{2k} \frac{R^2 - R_1^2}{\ln{(R/R_1)}} p_i. \end{equation}

At the outlet, the ambient pressure is set to zero, i.e.

(2.31)\begin{equation} p(x,y,z) = 0 \quad \text{for}\ \sqrt{x^2 + y^2} = R \quad \text{and} \quad 0 \leq z \leq h. \end{equation}

Furthermore, we exploit the rotational symmetry of the problem to save computational effort by simulating only a section of the domain, as shown in figure 2(b). This is done by applying the wedge boundary condition (see Greenshields & Weller (2022) for more details) at the wedge sides of the considered computational domain (cf. figure 3).

Figure 3. Discretized computational domain of the simplified simulation model. The image is scaled in the $z$ direction to improve the visibility.

2.2.3 Simulation set-up

The presented problem is solved numerically, using the unstructured finite volume method within OpenFOAM, applying the simpleFoam solver for steady-state single-phase Navier–Stokes equations that utilizes the SIMPLE algorithm developed by Caretto et al. (1973). The input data of the simulation set-up are publicly available, as well as the post-processing scripts and the secondary data (Pelz et al. 2023; Asghar et al. 2024).

Due to rotational symmetry of the problem, we consider only a slice of the whole three-dimensional domain (figure 2b). The slice is discretized using a prismatic (wedge) discretization. The angle of the wedge is $1^{\circ }$ in this study, and the depth of the domain running along the plane of symmetry comprises a single cell, as shown in figure 3. Parameters $N_r$ and $N_z$ denote the number of cells in radial and axial directions, respectively.

3. Results and discussion

In this section, a convergence study is performed to evaluate the accuracy of the simulation model. Below, the results of both modelling approaches are presented and compared with the measured data. The parameters used for the simulation studies are presented in table 1. The slip length value $\lambda = \lambda _o = \lambda _u = 540\ \textrm {nm}$ results from fitting experimental measurements to the analytical model (1.5) (see Corneli (2022) for details). In the following, we aim to validate the analytical model with numerical simulations. Since the inlet is composed of the same material, we will assume that $\lambda _i = \lambda = 540\ \textrm {nm}$.

Table 1. Parameters of the simulation model.

3.1 Mesh convergence study

In the following study, the torque on the bottom plate $M$ is calculated numerically for different mesh sizes, presented in table 2. The mesh width $\varDelta$ is defined by averaging the mesh widths $\varDelta _r$ and $\varDelta _z$ in the radial and axial directions, respectively, and is given as

(3.1)\begin{equation} \varDelta = \sqrt{\varDelta_r^2 + \varDelta_z^2} = \sqrt{\left(\frac{R}{N_r}\right)^2 + \left(\frac{h}{N_z}\right)^2}. \end{equation}

Furthermore, the error $e$ in the torque calculations is defined as

(3.2)\begin{equation} e = |M-M_{{ref}}|, \end{equation}

where the reference torque $M_{{ref}}$ is the torque resulting from the simulation with the finest-resolution mesh M6 (see table 2).

Table 2. Parameters used in the mesh convergence analysis.

The results of the convergence study are illustrated in figure 4. As expected, the results demonstrate a second-order convergence rate for the error $e$, given that a second-order discretization method is employed to solve the Navier–Stokes equation. Additionally, the torque $M$ is computed as the integral of the local velocity gradient over the area of the lower plate (see (2.23)).

Figure 4. Convergence analysis for the (simplified) simulation model.

3.2 Simulation results and validation

This section compares the velocity profiles obtained through simulation models with those derived from the analytical model. Moreover, it includes a comparative analysis of the torque values computed from the analytical model, simulations and experimental studies. This section aims to provide insights into the reliability of the models introduced in § 2. The results presented in this section are obtained using mesh M6 (see table 2) and a tribometer plate gap height $h=5\ \mathrm {\mu }$m. The values from the analytical model are considered as reference values.

3.2.1 Velocity profiles

Figure 5 shows the relative error, $e = (v-v_{{ref}})/v_{{ref}}$, in circumferential velocity $v$ along the radial direction for three different values of $z$. Both simulation models produce similar results, exhibiting that the inlet tube has no impact on the circumferential velocity within the domain. However, the analytical and simulation models differ in the inflow area, as depicted in figure 6. This deviation can be expected because the inflow is not modelled in the analytical model. Moreover, the deviation does not affect the torque at the bottom plate, since the torque depends only on the partial derivative ${\partial v}/{\partial z} |_{z=0}$.

Figure 5. Relative error, $e = (v-v_{{ref}})/v_{{ref}}$, in the circumferential velocity $v$ along the radial direction for actual and simplified simulation models.

Figure 6. Comparison of the circumferential velocity $v$ along the radial direction in the inflow region. Even with the presence of an inlet pipe with vanishing axial velocity, the flow almost immediately recovers the analytical model for the circumferential velocity after exiting the inlet.

Figure 7 shows the relative error $e$ in the velocity $v$ calculated along the axial coordinate $z$ for three different values of $r$. While the simulation results are in excellent agreement with the analytical model, a difference between the two simulation models can be observed. Additionally, the root-mean-square (RMS) error in circumferential velocity is presented in table 3.

Figure 7. Relative error $e$ in the circumferential velocity $v$ along the axial direction for actual and simplified simulation models.

Table 3. The RMS error in circumferential velocity $v$ along the radial and axial directions for simulation models with respect to the analytical model.

3.2.2 Torque

The available experimental measurements of torque $M$ and the tribometer plate gap height $h$, along with the operating parameters listed in table 4, are used to validate the models for torque measurement. The relationship of $M^{-1}$ and $h$ can be very well represented by the linear relation $M^{-1}=bh+a$. The parameters $a$ and $b$ are determined by Corneli (2022) using the least-squares method.

Table 4. Operating parameters of the SLT used to validate the models.

The analytical model measures the slip lengths of the upper and lower tribometer plates, $\lambda _u$ and $\lambda _o$, respectively, and the dynamic viscosity $\mu$ (see (2.23)). Subsequently, the following applies:

(3.3)\begin{equation} \left.\begin{array}{c} \lambda_u = \lambda_o = \dfrac{a}{2b}, \\ \mu = \dfrac{1}{b \varOmega I_p}. \end{array}\right\} \end{equation}

Furthermore, these physical parameters are then used to evaluate the analytical and simulation models.

The measurements of the tribometer plate gap height $h$ and the moment $M$ are subject to errors. In addition, temperature changes in the lubrication gap or variations in motor frequency can distort the measurements. The error propagation of the mentioned errors results in the uncertainties $\delta a$ and $\delta b$ of the parameters $a$ and $b$. From this, the $95\,\%$ confidence interval is derived, within which the measuring points are located with a probability of $95\,\%$.

Figure 8 shows the $95\,\%$ confidence interval of the measured data and the corresponding results of the analytical and simulation models. It is noticeable that the results of both modelling approaches agree very well with the measurement data, showing that the use of the analytical model is justified for evaluating the slip length from the measurements.

Figure 8. Validation of the models with the measurement of the torque $M$ as a function of the gap height $h$.

The effect of the velocity component $v_r$ in the radial direction on the torque $M$ is also investigated. Figure 9 shows the difference in $v_r$ obtained from actual and simplified simulation models at the outlet for an increase of the inlet pressure $p_i$. The behaviour of $v_r$ and $v_\phi$ shows that the simplified simulation model is a good approximation of the actual model. Figure 10 shows that the calculated torque is independent of the inlet flow pressure and, consequently, the radial velocity component.

Figure 9. Error in the radial velocity component at the outlet for an increase in the inlet pressure $p_i$, obtained from the actual and simplified simulation models.

Figure 10. Difference of reciprocal of the calculated torque $M$, using the actual and simplified models, with varying kinematic pressure at the inlet for gap height $h=5\ \mathrm {\mu }$m.

In this section, we verified the simplified analytical approach, which neglects the inertial effects as well as the feed flow entering and leaving the device during measurement of the torque. We show that the simplified model is sufficient to infer the slip length $\lambda$ from the experimentally measured torque.

4. Summary

This work verified and validated the simplified analytical model underlying the measurement principle of the SLT. The SLT is a classical plate tribometer developed to measure the effect of surface roughness in the form of slip length of the Navier-slip boundary condition. The simulation results show very good agreement between the experimental measurements, the analytical model and the simulation models, justifying the use of the simplified analytical model as the basis of the measurement principle of the SLT, which does not take into account the inlet and outlet boundary conditions as well as inertial effects.