# Optimization of mesh density for numerical simulations of hydrodynamic lubrication considering textured surfaces

Mohand Adjemout, Noël Brunetière, Jean Bouyer

► To cite this version:

Mohand Adjemout, Noël Brunetière, Jean Bouyer. Optimization of mesh density for numerical sim- ulations of hydrodynamic lubrication considering textured surfaces. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, 2015, 229 (9), pp.1132-1144. 10.1177/1350650115574535 . hal-04964015

HAL Id: hal-04964015

https://hal.science/hal-04964015v1

Submitted on 16 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimization of mesh density for numerical simulations of hydrodynamic lubrication considering textured surfaces

M. ADJEMOUT, N. BRUNETIERE, J. BOUYER

Institut Pprime, CNRS - Université de Poitiers - ENSMA, UPR 3346, Département Génie mécanique et systèmes Complexes, SP2MI, 11 Boulevard Marie et Pierre Curie, BP 30179, F86962 Futuroscope Chasseneuil Cedex, France

## Abstract

In this work, a parametric study of a lubricated textured contact is carried out using a numerical model. The model solves the Reynolds equation coupled with a mass-conserving cavitation algorithm, and addresses the hydrodynamic lubrication of a mechanical seal with textured surfaces. Three shapes of texture (square, spherical cap and triangle) are used, and the effect of the aspect ratio is also discussed for the case involving a triangular dimple. In order to estimate quantitatively the discretization error, the GCI (Grid Convergence Index) method is applied, considering several operating parameters (leakage, cavitation, load carrying capacity) and using three grid levels: coarse (100x100), medium (200x200) and fine (400x400). As could be expected, GCI analysis showed that there was a reduction in the discretization error when the grid system was refined. However, the study also showed that medium grid (200x200) solution has a GCI of less than 6%. Comparison between fine and medium mesh densities, with regard to accuracy and computational time requirements, indicated that a medium mesh density was appropriate in the case of the present program, since the computational time is reduced by a factor of nine when compared with the highest mesh density. Moreover, the size of the optimal mesh also depends on the operating conditions. The higher the lubrication number, the higher the number of nodes.

## 1 Introduction

Improving the tribological performance of lubricated contacts is a priority in relation to the reliability of mechanical components. Surface texturing is one of the main methods proposed and discussed in the literature to achieve this goal. It has been experimentally proven that surface texturing can significantly improve the performance of hydrodynamic lubricated contacts such as mechanical seals, bearings, etc. [1]. Surface texturing consists in modifying the topography of the mating surfaces by machining texture
micropatterns. The positive effect of textures is due to the hydrodynamic pressure buildup at the leading edge of texture patterns. This is associated with the appearance of cavitation at the trailing edge [2]. In order to understand the effect of surface texturing, several lubrication theories and approaches were developed and presented [3, 4, 5]. From most of the work to be found on the topic of lubricated contact, area density, depth ratio, shape and distribution of textures are presented as the important geometrical parameters which most influence the tribological performance. It is well known that textured surfaces generate several discontinuous zones on the mating faces. These discontinuous zones should be treated by the meshing. Thus, the usual solution adopted to solve this problem is to employ variable and adaptive meshes [6]. However, this method is difficult to implement due to the large number of texture patterns on the mating faces. Another method consists in using a regular fine grid to minimize the errors, but this solution has an impact on the computational time. Etsion et al [4] investigated numerically and experimentally the effect of spherical dimples on the load carrying capacity of mechanical seals. The authors noted that the best accuracy of their numerical simulations is obtained when the grid applied on the pores is five times greater than that applied outside of the pores. Arghir et al [7] analyzed theoretically the laminar flow between parallel surfaces, one of the surfaces having a regular macro-roughness pattern. They took a special care to eliminate any dependency of their numerical solutions on the grid by employing GCI (grid convergence index) for reporting grid independency. The same method (GCI) was employed by Dobrica et al [8] to analyze grid independency of a 3D THD (ThermoHydrodynamic) model with textured surfaces. Cupillard et al [9] applied the Richardson extrapolation to carry out an estimation of the error induced by the discretization spacing of a CFD simulation with a finite textured bearing. The same authors realized a similar study for thermal influence on textured inlet contact [10]. Qiu et al [11] theoretically studied the effect of various texture shapes on the load carrying capacity of a parallel slider bearing lubricated with gas. Selecting a Cartesian grid with 251 x 251 nodes per unit cell, the authors stated that the errors were inferior to 2%.

As a general comment, it can be said that all discrete methods introduce discretization error in the solution which is due to the choice of the space resolution (level of mesh). Despite the work cited previously [7-11], Estimating discretization error is not often considered in computational lubricated contact community even though surface texture is a hot topic in tribology [1]. Contrary to this, in the computational fluid dynamic community, discretization error estimation is mandatory to be allowed to publish numerical results in many
CFD Journals, as the case with ASME Fluids Engineering journal [12]. Despite the fact that numerical error cannot be eliminated completely in a simulation, it can be minimized or limited. For this purpose, an evaluation of the discretization errors is necessary to optimize the numerical program as well as the total number of nodes.

In this paper, the error in a numerical model of hydrodynamic lubrication for textured mechanical seals is evaluated. In this study, the authors focus on the discretization error. In order to estimate quantitatively the error resulting from the discretization, the GCI (Grid Convergence Index) method is used. This method is usually applied in CFD (Computational Fluid Dynamics) computation and is based on the generalized Richardson extrapolation theory [13]. This method has been recently recommended for finite element applications in solid and structural mechanics [14]. It enables an estimation of the discretization error to be made in the absence of knowledge of the analytical solution to the problem. In this study, the performance of a textured mechanical seal (load carrying capacity, leakage and cavitation ratio) has been computed and is presented in a dimensionless form.

## 2 Numerical model

The numerical model used in this study is based on the Reynolds equation coupled with a mass-conservative algorithm which deals with the film breakdown and reformation, generally known as cavitation. The cavitation model, based on the JFO model (Jacobsson-Floberg-Olsson), is derived from those published by Elrod [15], Payvar and Salant [16].

The geometrical configuration of the model is shown in (Figure 1(a)). Assuming that the lower ring is fixed, the Reynolds equation for a Newtonian fluid in a laminar flow regime is given by:

$$F \frac{\partial}{\partial r} \left( \frac{rh^3}{\mu} \frac{\partial D}{\partial r} \right) + F \frac{\partial}{\partial \theta} \left( \frac{h^3}{\mu r} \frac{\partial D}{\partial \theta} \right) = 6V \left[ \frac{\partial h}{\partial \theta} + (1-F) \frac{\partial}{\partial \theta} (hD) \right] \quad (1)$$

where F and D are the variables relating to the cavitation algorithm.
$$ \begin{cases} D = p - p_{cav} > 0 \\ \rho = \rho_0 \\ F = 1 \end{cases} \quad \begin{cases} p = p_{cav} \\ D = \frac{\rho}{\rho_0} - 1 < 0 \\ F = 0 \end{cases} \qquad (2) $$

$F$ is a switch variable indicating if cavitation occurs or not: for full film condition, $F$ equals 1, $D$ thus represents pressure, and equation (1) is Reynolds equation; When cavitation occurs, $F$ equals 0, $D$ represents the density and equation (1) becomes the mass conservation equation. The problem is supposed to be stationary, because the upper ring is perfectly smooth while the lower ring is textured. The studied domain is an angular portion of the ring with $R_i$ the inner radius and $R_0$ the outer radius, as shown in Figure1 (a). The finite volumes method is used to discretize the domain. For a control volume (Figure 2), the continuity equation is as follows:

$$ q_w - q_e + q_s - q_n = 0 \qquad (3) $$

This equation represents the balance of the incoming and outgoing mass flow rates of the control volume. The flow rates at the west and east sides can be expressed by:

$$ q_w = -\frac{h_w^3}{12\mu} R_w \Delta\theta \frac{\partial}{\partial r} (FD)_w \approx -\frac{h_w^3}{12\mu} R_w \Delta\theta \frac{F_p D_p - F_w D_w}{\Delta R} \qquad (4) $$

$$ q_e = -\frac{h_e^3}{12\mu} R_e \Delta\theta \frac{\partial}{\partial r} (FD)_e \approx -\frac{h_e^3}{12\mu} R_e \Delta\theta \frac{F_E D_E - F_p D_p}{\Delta R} \qquad (5) $$

At the south and north sides the flow rates are:

$$ q_s = -\frac{h_s^3}{12\mu} \frac{\Delta R}{R_s} \frac{\partial}{\partial \theta} (FD)_s + 0.5h_s \Delta R R_s \omega [1 - D_s(1 - F_s)] \approx -\frac{h_s^3}{12\mu} \frac{\Delta R}{R_s} \frac{F_p D_p - F_N D_N}{\Delta \theta} + 0.5h_s \Delta R R_s \omega [1 - D_s(1 - F_s)] \qquad (6) $$

$$ q_n = -\frac{h_n^3}{12\mu} \frac{\Delta R}{R_n \Delta \theta} \frac{\partial}{\partial \theta} (FD)_n + 0.5h_n \Delta R R_n \omega [1 - D_p(1 - F_p)] \approx -\frac{h_n^3}{12\mu} \frac{\Delta R}{R_n} \frac{F_N D_N - F_p D_p}{\Delta \theta} + 0.5h_n \Delta R R_n \omega [1 - D_p(1 - F_p)] \qquad (7) $$

More details on this approach can be found in references [17] and [18]. Equation (3) is solved using a direct LU decomposition method, taking into account the boundary conditions at the inner and outer radii, whereas a periodic condition is applied in the circumferential direction to ensure the continuity of the domain.
Finally, the simulation results give the pressure and the density fields after stabilization of the cavitation zone.

<figure><img src="image_2.png" /><figcaption>Figure 1: Configuration of the problem</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Finite volume method [17]</figcaption></figure>

In order to validate the current model, a comparison was carried out with the results published by Zheng and Meng [19]. The authors compared numerical results and experimental results obtained with a grooved ring immersed in oil at atmospheric pressure, while applying a normal load $F_N$ of 44 N. The data used in the simulation are given in Table 1. The results obtained for the film thickness at equilibrium are displayed in Figure 4 as a function of the bearing number. This figure shows that the results computed with the current model are in good agreement with the numerical results obtained by Zheng and Meng [19]. When compared to the experimental results, the numerical results exhibit a slight shift. It has to be noted that the measurement of film thickness is sometimes marred by a degree of error.
<table><caption>Table1: Parameters used in reference [19]</caption><thead><tr><th>Parameters</th><th>R<sub>i</sub> [m]</th><th>R<sub>o</sub> [m]</th><th>p<sub>i</sub> [Pa]</th><th>p<sub>o</sub> [Pa]</th><th>p<sub>cav</sub> [Pa]</th><th>F<sub>N</sub> [N]</th><th>μ [Pa.s]</th><th>ω [rpm]</th><th>γ [m]</th><th>δ</th></tr></thead><tbody><tr><th>Values</th><td>0.024</td><td>0.032</td><td>1.0×10<sup>5</sup></td><td>1.0×10<sup>5</sup></td><td>0.3×10<sup>5</sup></td><td>44</td><td>0.083</td><td>71-296</td><td>4.6×10<sup>-6</sup></td><td>56%</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure3: Results from the present model compared to those of reference [19].</figcaption></figure>

## 3 Discretization error and uncertainty estimation

Numerical errors that can occur in a simulation code are round-off error, iterative error and discretization error [20]. Round-off error occurs due to the use of finite numbers on digital computers. For example, in digital computation the number “1.0” is stored as floating number (0.99999...). Thus, its precision depends on the number of digits used after the point to the storage in computer’s memory. To mitigate the adverse effects of round-off error a double precision computation is considered in this study. Iterative error occurs as a consequence of the iterative method being applied to solve the discretization. It is defined as the deference between the numerical solution at a given iteration and the exact solution to the discrete equations. Thus, iterative error arises due to incomplete iterative convergence of a discrete system, for more information regarding iterative error see [21]. In the current model, as previously mentioned, a direct method is used to solve the discretized equation. The last source of error is the discretization error, this being defined (assuming zero round-off and iterative error) as the difference between the solution of the discrete equation and the analytically exact solution to the problem (analytical solution).
$$ \mathcal{E}(H) = f(H) - f_{\text{exact}} \qquad (8) $$

Discretization error is strongly related to the spatial resolution of the mesh. Thus, due to the limitation in computational resources, the discretization error is generally the largest error and the most difficult to estimate and eliminate [22]. When the exact solution is known, the traditional error analysis consists in finding a relation between the exact solution $f_{\text{exact}}$ and the numerical solution $f(H)$ as follows:

$$ f(H) = f_{\text{exact}} + C H^n + O(H^{n+1}) \qquad (9) $$

where $H$ is the cell size of the mesh, $C$ is a constant, $n$ is the order of accuracy (order of convergence) and $O(H^{n+1})$ is the higher order term involving powers of $H$. For $H \approx 0$, the higher order term $O(H^{n+1})$ can be neglected and the discretization error can be approximated by:

$$ \mathcal{E}(H) = f(H) - f_{\text{exact}} \approx C H^n \qquad (10) $$

Which, in logarithmic form, gives:

$$ \log(\mathcal{E}(H)) = \log(C) + n \log(H) \qquad (11) $$

However, the exact (analytical) solution is not available, which makes the discretization error more difficult to estimate. One of the traditional methods commonly used to estimate discretization error is the Richardson extrapolation. This method was developed in 1910 by L.F. Richardson [23] and consists to estimate extrapolated solution by combining solution on fine grid spacing ($H_1$) and coarse grid spacing ($H_2$) with $H_2/H_1=2$. Assuming a second numerical scheme, Richardson states that the solutions can be written as:

$$ \begin{cases} f_1 = f_{\text{extra}} + C H_1^2 + O(H_1^3) \\ f_2 = f_{\text{extra}} + C H_2^2 + O(H_2^3) \end{cases} \qquad (12) $$

The idea is then to combine the two solutions to eliminate the leading order error and then the constant $C$ is solved. Substituting this in to equation (12) a Richardson extrapolation solution $f_{\text{extra}}$ is obtained. At the end the Richardson extrapolated solution is expressed:
$$f_{\text{extra}} = f_1 + \frac{f_1 - f_2}{3} \qquad (13)$$

In the literature, numerous methods have been developed to estimate the discretization error. Most of these were developed for the finite elements method [24], but there are also more general methods of estimating the discretization error. These methods can be classified into two main categories. The first offers an estimate of the exact solution of the differential equation and includes the order refinement [25], the Richardson extrapolation [23, 26] and a recovery method from finite elements [27]. The second category of error estimator is based on the residual error and includes an error transport equation [28], finite residual methods [29] and adjoint methods [30]. A good summary of error estimators was given by C. J. Roy [31].

In the present work the authors have chosen the Grid Convergence Index method (GCI). Indeed, based on the generalized Richardson extrapolation [13], the GCI can be applied in any discretization scheme (finite elements, finite differences and finite volumes) and is relatively easy to implement without requiring an additional code. The GCI represents the absolute value of the Richardson extrapolation error multiplied by a factor. This method was firstly proposed by Roache [13] for uniform grids in computational fluid dynamics (CFD). Since then, several uncertainty estimators have been developed: the average local order method [32], the correction factors method [33], the factor safety method [34] and the least squares method [35]. Most of these are based on the GCI formula and are considered to be variants of the GCI method.

## 4 Grid Convergence Index (GCI)

This method requires multiple levels of mesh refinement to compute the extrapolated solution ($f_{\text{extra}}$) corresponding to zero cell size of mesh (this being an estimate of the exact solution). The value of the GCI indicates how far the solution is from this asymptotic solution. To apply this method the authors have followed the approach suggested by Celik et al [12].

Using equation (9) and replacing $f_{\text{exact}}$ by $f_{\text{extra}}$, the three discrete solutions for the cell size values $H_1, H_2$ and $H_3$, (with $H_1<H_2<H_3$) can be written for the three results of simulation, respectively $f_1, f_2$ and $f_3$,
$$
\begin{equation}
\begin{cases}
f_1 = f_{\text{extra}} + C H_1^n + O(H_1^{n+1}) \\
f_2 = f_{\text{extra}} + C H_2^n + O(H_2^{n+1}) \\
f_3 = f_{\text{extra}} + C H_3^n + O(H_3^{n+1})
\end{cases}
\tag{14}
\end{equation}
$$

If all calculations are in the asymptotic range, the higher order terms can be neglected. Solving the system of
equations (14), an order of accuracy can be computed:

$$
n = \frac{\left| \ln \frac{\varepsilon_{3,2}}{\varepsilon_{2,1}} \right| + \Theta(n)}{\ln(r_{2,1})} \quad (15)
$$

where

$$
\Theta(n) = \ln \left( \frac{r_{2,1}^n - z}{r_{3,2}^n - z} \right) \quad \text{with} \quad z = \operatorname{sign} \left( \frac{\varepsilon_{3,2}}{\varepsilon_{2,1}} \right) \quad \text{and where} \quad \varepsilon_{3,2} = f_2 - f_3, \varepsilon_{2,1} = f_1 - f_2 \tag{16}
$$

$r_{2,1}=H_2/H_1$ and $r_{3,2}=H_3/H_2$ are the refinement factors. In the case of a two dimensional model, these factors
can be expressed as a function of the total number of nodes N (with $N_1>N_2>N_3$):

$$
r_{2,1} = \sqrt{\frac{N_1}{N_2}} \tag{17a}
$$

$$
r_{3,2} = \sqrt{\frac{N_2}{N_3}} \tag{17b}
$$

Equation (15) is a transcendental equation and can be solved with an iterative method. However, for a
constant refinement factor r<sub>2,1</sub>=r<sub>3,2</sub>, the Θ(n) term is equal to zero, and equation (15) becomes:

$$
n = \frac{\left| \ln \frac{\varepsilon_{3,2}}{\varepsilon_{2,1}} \right|}{\ln(r_{2,1})} \quad (18)
$$

If the order of accuracy *n* is higher than the theoretical order of the method, this may correspond to a super-
convergence [35]. Therefore, in this case, *n* is replaced by the theoretical order of the numerical scheme. The
theoretical order of the numerical scheme for the present program is two (centered scheme). When the order
of accuracy is obtained, an estimate of the relative extrapolated errors for the fine and medium grids can be computed using the expressions:

$$ \frac{f_{extra} - f_1}{f_1} \approx \frac{e_{2,1}}{r_{2,1}^n - 1} \qquad (19a) $$

$$ \frac{f_{extra} - f_2}{f_2} \approx \frac{e_{3,2}}{r_{3,2}^n - 1} \qquad (19b) $$

where $e_{2,1}$ and $e_{3,2}$ are the relative errors

$$ e_{2,1} = \left| \frac{f_1 - f_2}{f_1} \right| \qquad (20a) $$

$$ e_{3,2} = \left| \frac{f_2 - f_3}{f_2} \right| \qquad (20b) $$

The extrapolated solution $f_{extra}$ can be derived from equations (19a) or (19b). However, an accurate estimation of the extrapolated solution corresponds to the simulations with fine grids spacing i.e with equation (19a) where $H_1.<H_2$ or ($N_1>N_2$).

As previously explained, the GCI represents the absolute value of the Richardson extrapolation multiplied by a factor of safety. Thus, the GCI at fine grid and medium grid can be expressed by:

$$ GCI_{2,1} = F_s \frac{e_{2,1}}{r_{2,1}^n - 1} \qquad (21a) $$

$$ GCI_{3,2} = F_s \frac{e_{3,2}}{r_{3,2}^n - 1} \qquad (21b) $$

$F_s$ is a safety factor which ensures increased confidence in the estimation of uncertainty. According to Roache [36], when three discrete solutions are available, the safety factor is set at 1.25. This value is also used by Celik et al [12]. When the values of $GCI_{3,2}$ and $GCI_{2,1}$ are obtained, it is then possible to check the condition:
$$ \frac{GCI_{3,2}}{GCI_{2,1}} \approx r_{2,1}^n \qquad (22) $$

If this condition is satisfied, it indicates that the discrete solutions are in the asymptotic range. It should be noted that the Richardson extrapolation is reliable only for solutions respecting the monotonic convergence condition given by $0 < \varepsilon_{2,1}/\varepsilon_{3,2} < 1$. This condition was generalized in the work of Celik et al [37] for the case of oscillatory convergence. By introducing the absolute value in the order of accuracy previously shown in equation (15) and (18), the condition becomes

$$ 0 < |\varepsilon_{2,1}/\varepsilon_{3,2}| < 1 \qquad (23) $$

This method cannot be used when $\varepsilon_{3,2}$ or $\varepsilon_{2,1}$ are very close to zero ($10^{-15}$ in this work, double precision being used in the computations). This may either indicate that the exact solution has been reached or give an indication of oscillatory convergence. As a consequence, throughout this study, all similar cases are excluded from the results of GCI. They generally correspond to cases where no cavitation occur leading to high convergence even with a coarse mesh.

## 5 Grid convergence study

Three shapes of dimples are used on the stationary ring with the orientation shown in Figure 4. Each dimple is entrapped in a “square cell” with a length of L=200 µm. It should be noted that, since the cell is a part of a disk, its shape is nearly a “square” with a weak curvature since the ratio of its width to the radius is $\Delta R/R_i=0.0077$. Referring to the data given in table 2, each ring contains 808 dimples on its circumference.

<figure><img src="image_5.png" /><figcaption>Figure4: Textured ring (exaggerated size)</figcaption></figure>
Table 2 summarizes the data used in the simulation, while Table.3 presents the cases considered in this study.

<table><caption>Table 2: Data used in the simulations</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>R<sub>i</sub> [m]</td><td>2.575×10<sup>-2</sup></td></tr><tr><td>R<sub>o</sub> [m]</td><td>2.595×10<sup>-2</sup></td></tr><tr><td>p<sub>i</sub> [Pa]</td><td>0.10×10<sup>5</sup></td></tr><tr><td>p<sub>o</sub> [Pa]</td><td>0.10×10<sup>5</sup></td></tr><tr><td>p<sub>cav</sub> [Pa]</td><td>0</td></tr><tr><td>μ [Pa.s]</td><td>1.0×10<sup>-3</sup></td></tr><tr><td>ω [rpm]</td><td>25-250</td></tr><tr><td>δ</td><td>30%</td></tr><tr><td>L [m]</td><td>2.0×10<sup>-4</sup></td></tr><tr><td>h<sub>f</sub> [m]</td><td>1.0×10<sup>-6</sup></td></tr></tbody></table>

<table><caption>Table 3: Cases tested in this study</caption><thead><tr><th>Nodes</th><th colspan="4">Texture patterns</th></tr></thead><tbody><tr><td>100X100</td><td rowspan="3">L<br>ΔV<br>L</td><td rowspan="3">L<br>ΔV<br>L</td><td rowspan="3">L<br>ΔV<br>L</td><td rowspan="3">L<br>ΔV<br>L</td></tr><tr><td>200X200</td></tr><tr><td>400X400</td></tr><tr><td>factor λ</td><td>1</td><td>5/2</td><td>5</td><td>5/2</td><td>5/2</td></tr></tbody></table>

In order to analyze the influence of the discretization error on the simulation results, three performance parameters are analyzed:

### Dimensionless load carrying capacity

$$ \bar{p} = \frac{1}{S \times p_0} \int_S p \cdot ds - 1 \qquad (24) $$

### Average radial leakage

$$ \bar{Q} = \frac{12\mu(R_0 - R_i)}{2\pi p_0 h_f^3 R_m} Q \qquad (25) $$
Cavitation Ratio

$$ \bar{\epsilon} = \frac{N_c}{\delta \times N_i} \qquad (26) $$

All the three parameters will be presented as a function of a dimensionless hydrodynamic parameter which is given by:

$$ \psi = \frac{6\mu\omega R_m (R_0 - R_i)}{p_o h_f^2} \qquad (27) $$

6 Results and discussion

Figure 5 shows examples of pressure distributions for three dimple shapes for $\lambda=5/2$ and $\psi = 776$. The pressure contour shows the cavitation occurrence predicted by the JFO model. It reveals the two mechanisms of cavitation: the film breakdown at the upstream rim of the dimple and the film reformation at the downstream rim of the dimple. The film reformation is followed by a pressure build up in the sliding direction. It is known that the pressure field has a dependency on the operating conditions and on the type of lubricant. In addition, Figure 5 shows that the size of the cavitation zone and the magnitude of the pressure peak obviously depend on the dimple shape.

<figure><img src="image_6.png" /><figcaption>Figure 5: Single dimple pressure distribution ($\lambda=5/2$, $\psi=776$)</figcaption></figure>

During this analysis, the results of numerical uncertainties are presented as error bars using the GCI with an average value of the order of accuracy $n = n_{av}$, as suggested by Celik et al [12]. The maximum values of GCI observed over the range of the hydrodynamic parameter are displayed in Table.4.
<table><caption>Table 4: Maximum values of GCI obtained for each case ($N_1=400 \times 400$, $N_2=200 \times 200$, $N_3=100 \times 100$)</caption><thead><tr><th>Shape</th><th colspan="3">Square</th><th colspan="3">Spherical cap</th><th colspan="6">Triangle</th></tr><tr><th>Ratio λ</th><th colspan="3">5/2</th><th colspan="3">5/2</th><th colspan="3">1</th><th colspan="3">5/2</th><th colspan="3">5</th></tr><tr><th>variable</th><th>p&#773;</th><th>Q&#773;</th><th>ζ&#773;</th><th>p&#773;</th><th>Q&#773;</th><th>ζ&#773;</th><th>p&#773;</th><th>Q&#773;</th><th>ζ&#773;</th><th>p&#773;</th><th>Q&#773;</th><th>ζ&#773;</th><th>p&#773;</th><th>Q&#773;</th><th>ζ&#773;</th></tr></thead><tbody><tr><td>n<sub>av</sub></td><td>1.86</td><td>0.63</td><td>1.15</td><td>0.61</td><td>1.00</td><td>1.00</td><td>0.93</td><td>0.93</td><td>1.00</td><td>1.23</td><td>0.62</td><td>1.40</td><td>1.95</td><td>0.81</td><td>0.83</td></tr><tr><td>GCI<sub>3,2</sub> %</td><td>0.29</td><td>88.9</td><td>11.1</td><td>0.51</td><td>92</td><td>2.3</td><td>1.87</td><td>9.23</td><td>3.05</td><td>1.28</td><td>5.91</td><td>6.86</td><td>0.33</td><td>5.51</td><td>14.8</td></tr><tr><td>GCI<sub>2,1</sub> %</td><td>0.1</td><td>76.4</td><td>2.53</td><td>0.36</td><td>72.4</td><td>1.07</td><td>0.45</td><td>7.56</td><td>1.31</td><td>0.53</td><td>2.59</td><td>3.23</td><td>0.08</td><td>2.84</td><td>9.56</td></tr></tbody></table>

The intention in this analysis is to describe the results relating to numerical uncertainties, as well as to address the mechanical seal performance. Hence, a short discussion of mechanical seal performance is given for each case in the results.

## 6.1 Influence of the dimple shape

<figure><img src="image_7.png" /><figcaption>Figure 6: Dimensionless leakage for three dimple patterns</figcaption></figure>

Figure 6 gives the results of dimensionless leakage ($\bar{Q}$) corresponding to the three shapes of dimple having a ratio $\lambda=5/2$. In this figure, the vertical axis on the left hand side corresponds to the leakage associated with square and spherical cap dimples, while the one on the right hand side corresponds to the leakage stemming from the triangular dimple. The triangular dimples give a high leakage value, about a thousand times greater in absolute value than the leakage arising when other dimple shapes are used. This high value is due to the pumping effect which is a consequence of the shape and the orientation of the triangular dimples. The negative sign of the leakage means that the fluid flows from the outer radius to the inner radius of the
mechanical seal. Regarding the error bars, the discretization uncertainty is very high for the cases of the square and spherical cap dimples. Here, the maximums values of GCI are respectively $GCI_{2,1}^{\max} = 76.4\% - 72.4\%$ and $GCI_{3,2}^{\max} = 88.9\% - 92\%$. This may indicate that the numerical solutions are not in the asymptotic range and are significantly influenced by the discretization error. These high uncertainties of discretization can be explained by the low leakage generated when these dimple shapes are used. This makes the results highly sensitive to the discretization error. Indeed, the common geometrical characteristic of the circular and square shape is that these two dimples are symmetrical with respect to the sliding direction. Thus, the pressure field is also symmetrical. Therefore, it is expected that the flow rate in the radial direction must be close to zero (it would be zero for a perfect square cell) because the flow rate at the inner radius must be close to flow rate at the outer radius due to the symmetrical shape of the pressure field. The numerical results show that, with these two dimples, the leakage is very small but is not equal zero. The leakage values are thus merged with the error induced by the discretization and very sensitive to the mesh refinement, leading to high GCI values. For the leakage observed with the triangle, the values of GCI are $GCI_{2,1} < 3\%$ and $GCI_{3,2} < 6\%$.

<figure><img src="image_8.png" /><figcaption>Figure 7: Dimensionless load carrying capacity for three dimple patterns</figcaption></figure>

Figure 7 shows the results of the dimensionless load carrying capacity ($\bar{P}$) for the three shapes of dimples. For the three shapes tested here, the spherical cap dimple exhibits a greater ability to enhance the load carrying capacity, while the triangular dimple is unable to generate a positive dimensionless load. The
maximum discretization uncertainties with respect to this variable are low. Grid convergence study shows that the discretization uncertainties are less than 1% with fine grid (400x400) and less than 2% with the medium grid (200x200).

<figure><img src="image_9.png" /><figcaption>Figure 8: Cavitation ratio for the three dimples</figcaption></figure>

Figure 8 shows the results for the variation in cavitation ratio ($\bar{\zeta}$) in relation to the shape of the dimples. As expected, the cavitation ratio increases with the hydrodynamic parameter. However, it remains at roughly the same level for the three dimple shapes, especially for high values of the hydrodynamic parameter. The error bars show that the discretization uncertainties are low and do not exceed 5%, except in the case where GCI = 11.1%, obtained with the medium grid (200x200) and $\psi = 77.6$.

More results of GCI estimator with more refined mesh are given in Appendix (Figure13- Figure18), for the load carrying capacity and the cavitation ratio of the three dimple shapes with a depth ratio $\lambda = 5/2$.

## 6.2 Influence of the depth ratio of the triangular dimple

The dimensionless leakages for the triangular dimples when $\lambda$ is varied are shown in Figure 9. The highest dimensionless leakage, expressed as an absolute value, is observed with a depth ratio of 5/2. It seems that an optimal value of depth ratio exists between $\lambda=1$ and $\lambda=5$. The discretization uncertainties were decreased when the depth ratio was increased.
<figure><img src="image_10.png" /><figcaption>Figure 9: Dimensionless leakage for triangular dimples in relation to varying $\lambda$ ratio</figcaption></figure>

The highest maximum values of GCI are $GCI_{3,2}^{\max}=9.23\%$ and $GCI_{2,1}^{\max}=7.56\%$, these being observed at a depth ratio of $\lambda=1$ and $\psi=776$. However, for the other results, the GCI were observed to be less than 6% for the medium grid and less than 3% for the fine grid. It can thus be stated that a discretization error has a slight influence on the solutions under consideration. Moreover, it should be noted that one result relating to the hydrodynamic parameter $\psi = 77.6$ was deleted from the presented results because its relative error was very close to zero.

<figure><img src="image_11.png" /><figcaption>Figure 10: Dimensionless load carrying capacity in relation to varying $\lambda$ ratio</figcaption></figure>

Regarding the results shown in Figure 10, the dimensionless load carrying capacity was improved with a depth ratio value of 5. However, it is still unable to give a positive dimensionless load carrying capacity. This
amelioration can be explained by the depth of dimple. Here, $\lambda$= 5 corresponds to a deeper dimple, a consequence being that the mating faces are better supplied with lubricant thus limiting the cavitation development. As for discretization uncertainty, most of results are in the asymptotic range, where $GCI_{3,2} < 2\%$ and $GCI_{2,1} < 1\%$. An exception was observed for the case corresponding to a hydrodynamic parameter $\psi$ =155 and depth ratio $\lambda$ =1. In this case, the generated load carrying capacity is close to zero even if large pressure peaks are observed, due to large cavitation zones. A slight variation in pressure can thus significantly modify the load carrying capacity, which leads to create a high dependence of this parameter on mesh density. Note that two simulation results associated with the hydrodynamic parameter value of 77.6 are omitted, their relative error being close to zero.

<figure><img src="image_12.png" /><figcaption>Figure 11: Cavitation ratio for triangular dimples in relation to varying $\lambda$ ratio</figcaption></figure>

Figure 11 shows that the cavitation ratio decreases with increasing depth ratio $\lambda$, whereas the discretization uncertainty remains more or less constant. The highest maximum values of $GCI$ with respect to this variable are observed with the depth ratio $\lambda$=5, those values being respectively $GCI_{3,2}^{\max}$ =14.8% and $GCI_{2,1}^{\max}$ = 9.56%, which are observed at a hydrodynamic parameter of $\psi$ =155. Apart from these two results, most of the results show a weak dependence of the cavitation on the grid refinement, where $GCI_{3,2} < 5\%$ and $GCI_{2,1} < 2\%$. Note again that one case is deleted from these results: that corresponding to the low value of the hydrodynamic parameter and a depth ratio $\lambda$= 5. In this case, the relative error is close to zero.

## 6.3 Computation time
In all the cases presented, the values of $GCI_{2,1}$ were lower than $GCI_{3,2}$. This means that the discretization error was logically minimized by reducing the size of cells. However, the computational time when reducing the cell size can become very high. For example, the discretization uncertainty of load carrying capacity and the corresponding computational times are shown in Figure 12. This figure relates to four levels of discretization and a hydrodynamic parameter $\psi=620$. The results of simulation with 50X50 nodes were added to afford an estimate of the discretization uncertainty when using 100X100 nodes. Examining the results presented in Figure 12, the computational time (see Figure12 (b)) for a set-up using 400X400 nodes is about nine times greater than for the 200X200 node set-up. By contrast, the discretization uncertainty decreases. Accordingly, the choice of mesh density is defined by a trade-off between the goal of accuracy and the

<figure><img src="image_13.png" /><figcaption>Figure 12: Discretization uncertainty and computational time versus mesh density:
(a) Discretization uncertainty, (b) Computational Time
(Load carrying capacity of spherical cap dimple at $\psi=6.20 \ 10^{+2}$)</figcaption></figure>
negative aspects of increased computational time.

## 7 Conclusion

In this study, an analysis of the influence of mesh density on the convergence of a numerical model of lubricated textured surfaces was carried out. Using the GCI method, it was found that the discretization error varies as a function of the dimple shape and the depth ratio $\lambda$ of the texture. It was also found that the discretization error depends on the selection of the studied (/calculated) quantity (load carrying capacity, cavitation and leakage). The results showed that the radial leakage was severely marred by the discretization error, especially with the square and spherical cap dimples. Apart from these results, the discretization error remains in an acceptable range for load carrying capacity and cavitation ratio. In the simulations carried out by the present authors, the highest mesh density of 400X400 gives the most accurate results: $GCI < 1\%$ for load carrying capacity, $GCI < 2\%$ for the cavitation ratio and $GCI < 3\%$ for the leakage (except the square and spherical cap). However, the computational time is longer with this mesh density. With a lower mesh density of 200X200, the results are slightly perturbed by the discretization error: $GCI < 2\%$ for the load carrying capacity, $GCI < 5\%$ for the cavitation ratio and $GCI < 6\%$ for the leakage. With this mesh density the computational time is reduced by a factor of nine when compared with the highest mesh density. In this study, only one dimple was considered in the studied domain (the radial band), while in practice, many patterns will have to be considered along the radial band. In order to keep the computational time and the accuracy of the results within an acceptable range, the density of 200X200 appears to be a good compromise between accuracy and computation time. Note that the Appendix of the paper offers a more detailed view of the GCI for a large spectrum of mesh density. It is thus possible to optimize the mesh for a target accuracy.
## Funding

The work presented in this study was supported by the French National Research Agency (grant no ANR-11-RMNP-0008)

## Nomenclature

<table><tr><td>r</td><td>= radial coordinate (m)</td></tr><tr><td>θ</td><td>= angular coordinate (rad)</td></tr><tr><td>h</td><td>= fluid film thickness (m)</td></tr><tr><td>h<sub>f</sub></td><td>= imposed film thickness (m)</td></tr><tr><td>μ</td><td>= dynamic viscosity (Pa.s)</td></tr><tr><td>V</td><td>= sliding velocity (m/s)</td></tr><tr><td>F</td><td>= switch variable</td></tr><tr><td>D</td><td>= universal variable</td></tr><tr><td>p</td><td>= local pressure (Pa)</td></tr><tr><td>p<sub>cav</sub></td><td>= cavitation pressure (Pa)</td></tr><tr><td>p<sub>i</sub>, p<sub>o</sub></td><td>= pressure at the inner and outer radii (Pa)</td></tr><tr><td>ρ</td><td>= local density (kg/m<sup>3</sup>)</td></tr><tr><td>ρ<sub>0</sub></td><td>= fluid density (kg/m<sup>3</sup>)</td></tr><tr><td>q</td><td>= local fluid mass flow (kg/s)</td></tr><tr><td>Q</td><td>= leakage rate (m<sup>3</sup>/s)</td></tr><tr><td>R<sub>i</sub>, R<sub>o</sub></td><td>= inner and outer radii of the seal (m)</td></tr><tr><td>R<sub>m</sub> = (R<sub>i</sub> + R<sub>o</sub>) / 2</td><td>= mean seal radius</td></tr><tr><td>ω</td><td>= rotating speed (rad/s)</td></tr><tr><td>S = π(R<sub>o</sub><sup>2</sup> - R<sub>i</sub><sup>2</sup>)</td><td>= sealing area (m<sup>2</sup>)</td></tr><tr><td>γ</td><td>= texture depth (m)</td></tr><tr><td>λ = h<sub>f</sub> / γ</td><td>= depth ratio</td></tr><tr><td>δ</td><td>= density of texture</td></tr><tr><td>F<sub>N</sub></td><td>= closing force (N)</td></tr><tr><td>W = F<sub>N</sub> / S</td><td>= Nominal mean contact pressure (Pa)</td></tr><tr><td>G = μω / W</td><td>= bearing number</td></tr></table>
$f_{exact}$ = exact solution

$f_{extra_{j+1,j}} = f_j + \frac{f_{j+1} \times e_{j+1,j}}{r_{j+1,j}^n - 1}$ = extrapolated solution

$N_j$ = total number of nodes

$N_c$ = number of cavitating nodes

$H_j$ = mesh cell size corresponding to the total number of nodes $N_i$

$n = \frac{|\ln \frac{\varepsilon_{3,2}}{\varepsilon_{2,1}}|}{\ln(r_{2,1})}$ = order of accuracy

$n_{av} = \frac{1}{Ns} \sum n$ = averaged order of accuracy

$Ns$ = number of simulations

$r_{j+1,j} = \sqrt{\frac{N_j}{N_{j+1}}}$ = refinement factor

$f$ = variable to be solved ($\bar{Q}, \bar{p}$ or $\bar{\zeta}$)

$f_j$ = numerical solution with mesh cell size $H_j$

$\varepsilon_{j+1,j} = f_j - f_{j+1}$ = difference between two solutions at two mesh size levels

$e_{j+1,j} = \left| \frac{f_j - f_{j+1}}{f_j} \right|$ = relative error

$GCI_{j+1,j} = F_s \frac{e_{j+1,j}}{r_{j+1,j}^n - 1}$ = grid convergence index

$GCI_{j+1,j}^{\max}$ = maximum values of $GCI_{j+1,j}$ observed over the range of $\psi$

$F_s$ = safety factor

$\bar{p} = \frac{1}{S \times p_0} \int_S P \cdot ds - 1$ = dimensionless load carrying capacity

$\bar{Q} = \frac{12 \mu (R_o - R_i)}{2 \pi \rho_0 h_f^3 R_m} Q$ = dimensionless radial leakage

$\bar{\zeta} = \frac{N_c}{\delta \times N_i}$ = cavitation ratio

$\psi = \frac{6 \mu \omega R_m (R_o - R_i)}{p_0 h_f^2}$ = hydrodynamic parameter
**Subscripts**

e,n,s,w = east (face), north, south, west

E,N,P,S = east (node), north, current, south, west

j = mesh level
References

1. Etsion I. State of the art in Laser surface texturing. *Transaction of the ASME* 2005; 127: 248-253.

2. Hamilton D B, Walowit J A and Allen C M. A theory of lubrication by micro-irregularities. *ASME J Basic Eng-T*, 1966; 88:177-185.

3. Etsion I and Burstein L. A model for mechanical seals with regular microsurface structure. *Tribol T* 1996; 39: 677-83.

4. Etsion I, Kligerman Y and Halperin G. Analytical and experimental investigation of laser-textured mechanical seal faces. *Tribol T* 1999; 42: 511-6.

5. Qiu Y and Khonsari M. On the prediction of cavitation in dimples using a mass- conservative algorithm. *ASME J Tribology*, 2009; 131: 041702-11.

6. Ali M S, Doolan C J and Wheatley V. Grid convergence study for a two dimensional simulation of flow around a square cylinder at a low Reynolds number. *Seventh International Conference on CFD in the Minerals and Process Industries CSIRO, Melbourne, Australia 9-11 December 2009*

7. Arghir M, Roucou N, Helene M et al. Theoretical Analysis of the Incompressible Laminar Flow in a Macro-Roughness Cell. *ASME J Tribology* 2003; 125: 309-318.

8. Dobrica M B, Fillon M. Thermohydrodynamic Behavior of a Slider Pocket Bearing. *ASME J Tribology* 2006; 128: 312-318.

9. Cupillard S, Cervantes M J and Glavatskih S. A CFD Study of a Finite Textured Journal Bearing. 24th Symposium on Hydraulic Machinery and Systems.

10. Cupillard S, Cervantes M J and Glavatskih S. Thermal Analysis of Lubricant Flow in a Textured Inlet Contact. *The 19th International Symposium on Transport Phenomena*, 17-20 August, 2008, ICELAND.

11. Qiu M, Delic A and Raeymaekers B. The Effect of Texture Shape on the Load-Carrying Capacity of Gas-Lubricated Parallel Slider Bearings. *Tribol Lett* 2012; 48:315-327.
12. Celik I B, Ghia U and Roche P J et al. Procedure for estimation and reporting of uncertainty due to discretization in CFD applications. *ASME J Fluid Eng* 2008; 130: 1-4.

13. Roache P J. Perspective: A method for uniform reporting of grid refinement studies. *ASME J Fluid Eng* 1994; 116:405-413.

14. Kwaśniewski L. Application of grid convergence index in FE computation. *Bulletin of the Polish Academy of Sciences, Technical Sciences* 2013; Vol. 61

15. Elrod H G. A cavitation algorithm. *J Lubr Technol* 1981; 103: 350-354

16. Payvar P and Salant R F. A computational method for cavitation in a wavy mechanical seal. *ASME J Tribology* 1992; 114:199-204.

17. Minet C, Brunetière N and Tournerie B. A deterministic mixed lubrication model for mechanical seals. *ASME J Tribology* 2011; 133: 042203-11.

18. Brunetière N and Tournerie B. Numerical analysis of a surface-textured mechanical seal operating in mixed lubrication regime. *Tribol int* 2012; 49: 80-89.

19. Zhang J and Meng Y. Direct observation of cavitation phenomenon and hydrodynamic lubrication analysis of textured surfaces. *Tribol Lett* 2012; 46: 147-158.

20. Roy C J. Review of code and solution verification procedures for computational. *Journal of Computational Physics* 205 (2005) 131-156.

21. Oberkampf W L and Roy C J. Verification and Validation in scientific computing. Cambridge University press, Cambridge 2010.

21. Oberkampf W L and Roy C J. Verification and Validation in scientific computing. *Cambridge University press, Cambridge 2010*

22. Phillips T S. Residual Methods for Discretization Error Estimation. Presented at 20th AIAA Computational Fluid Dynamics Conference 27-30, Honolulu, Hawaii, June 2011.
23. Richardson L F. The approximate arithmetical solution by finite differences of physical problems involving differential equations with an application to the stresses in a masonry dam. *Trans. Royal Society London* 1910 ; Ser. A 210 : 307-357.

24. Ainsworth M and Oden J T. A Posteriori Error Estimation in Finite Element Analysis. *Comput Method Appl M* 1997; 142: 1-88.

25. Fehlberg E. Low-order classical Runge-Kutta formulas with step size control and their application to some heat transfer problems, *NASA Technical Report 315, National Aeronautics and Space Administration* July 1969.

26. Richards S A, Completed Richardson extrapolation in space and time. *Commun Numer Meth En* 1997; 13: 573-582.

27. Zienkiewicz O C and Zhu J Z. The super convergent patch recovery and a posteriori error estimates, Part 2: error estimates and adaptivity. *Int J Numer Meth Eng* 1992; 33: 1365-1382.

28. Zhang X D, Trepanier J Y and Camarero R. A Posteriori error estimation for finite-volume solutions of hyperbolic conservation laws, *Comput Method Appl M* 2000; 185: 1-19.

29. Babuska I and Rheinboldt W C. Error estimates for adaptive finite element computations. *SIAM J Numer Anal* 1978; 15: 736-754.

30. Pierce N A and Giles M B. Adjoint recovery of super convergent functionals from PDE approximations, *SIAM Rev* 2000; 42: 247-264.

31. Roy C J. Review of discretization error estimators in scientific computing. *Presented at the 48th AIA Aerospace Sciences Meeting, 4-7, Orlando, Florida, January 2010.*

32. Cadafalch J, Perez-Segarra C D and R Consul et al. Verification of finite volume computations on steady-state fluid flow and heat transfer, *ASME J Fluid Eng* 2002; 124: 11-21.
33. Stern F, Wilson R and Shao J. Quantitative V&V of CFD simulations and certification of CFD codes. *Int J Numer Meth Fl* 2006; 50: 1335–1355.

34. Xing T and Stern F. Factors of safety for extrapolation. *IIHR Technical Report No. 469*, March 2009.

35. Eca L, Hoekstra M. Evaluation of numerical error estimation based on grid refinement studies with the method of manufactured solutions. *Comput Fluids* 38; 1580-1591.

36. Roache P J. Verification and Validation in computational science Engineering. *Hermosa Publishers, New Mexico*, 1998.

37. Celik I and Karatekin O. Numerical experiments on application of Richardson extrapolation with no uniform grids. *ASME J Fluid Eng* 1997; 119: 584-590.
**Appendix:** Extended discretization uncertainties (GCI) of the three dimple shapes for a depth ratio λ= 5/2

<figure><img src="image_14.png" /><figcaption>FIG.13- Discretization uncertainties on load of the triangle dimple</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG.14- Discretization uncertainties on the cavitation of the triangle dimple</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIG.15- Discretization uncertainties on load of the spherical cap dimple</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIG.16- Discretization uncertainties on the cavitation of the spherical cap dimple</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIG.17- Discretization uncertainties on the load of the square dimple</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIG.18- Discretization uncertainties on the cavitation of the square dimple</figcaption></figure>