A Local Adaptive Mesh Reﬁnement for JFO Cavitation Model on Cartesian Meshes

: Nonuniform mesh is beneﬁcial to reduce computational cost and improve the resolution of the interest area. In the paper, a cell-based adaptive mesh reﬁnement (AMR) method was developed for bearing cavitation simulation. The bearing mesh can be optimized by local reﬁnement and coarsening, allowing for a reasonable solution with special purpose. The AMR algorithm was constructed based on a quadtree data structure with a Z-order ﬁlling curve managing cells. The hybrids of interpolation schemes on hanging nodes were applied. A cell matching method was used to handle periodic boundary conditions. The difference schemes at the nonuniform mesh for the universal Reynolds equation were derived. Ausas’ cavitation algorithm was integrated into the AMR algorithm. The Richardson extrapolation method was employed as an a posteriori error estimation to guide the areas where they need to be reﬁned. The cases of a journal bearing and a thrust bearing were studied. The results showed that the AMR method provided nearly the same accuracy results compared with the uniform mesh, while the number of mesh was reduced to 50–60% of the number of the uniform mesh. The computational efﬁciency was effectively improved. The AMR method is suggested to be a potential tool for bearing cavitation simulation.


Introduction
Hydrodynamic bearings are widely used in machines, particularly in engines and power plants. The bearing consists of a rotating journal inserting a bore or a sleeve [1]. To ensure adequate lubrication, the journal wall driving velocity should be high enough to force oil into a converging clearance so that the load can be supported by the increasing hydrodynamic pressure [2]. The load-carrying capacity is one of the critical performance parameters and is treated as the prime objective in journal bearing design [1].
Cavitation can occur at particular kinds of bearings where the variable-shaped divergence clearance exists in the oil passage. Typically, cavitation occurs at the divergence clearance region in full circle bearings. Finger-shaped voids that cover most of the divergence clearance region can be observed [3]. Another specific case is textured thrust bearings (such as dimple-enhanced seal-like thrust bearings), in which cavitation occurs at the tail end of each dimple [4]. Although cavitation does not necessarily have a deleterious effect upon the load-carrying capacity, the cavitation algorithm's predicted load-carrying capacity is significantly affected [5].
Cavitation is essentially characterized as the two-phase flow of liquid and gas [6]. The liquid phase is oil. The gas phase is the air which escapes from oil when the pressure is below saturation pressure, or the oil vapor where the oil boils to form bubbles if the pressure is lower than the vapor pressure [7]. No matter what kind of cavitation, the simple data structures (permit reuse of uniform mesh code). However, it becomes challenging to encode as grid hierarchy increases. (2) Cell-based AMR refines only those cells that are supposed to be refined. The nature of structured meshes is lost due to independent cell management. However, a tree-based data structure is often used to manage the meshes, making the neighbor relation efficiently obtained [31]. The unstructured AMR also provides geometric flexibility at the cost of explicitly storing all neighborhood relations [32]. In the paper, a cell-based AMR algorithm was developed for bearing cavitation simulation. The aim is twofold: (1) The AMR method is introduced into cavitation problem. It provides another way to significantly improve the computational efficiency. This method exhibits potential advantages in analysis of bearings with complex geometric shapes. (2) Due to the complexity of the AMR method, the basic AMR algorithm was explained in the paper. The algorithm may provide guidance for people interested in AMR. The AMR algorithm was constructed based on a quadtree data structure with a Zorder filling curve managing cells. The hybrids of interpolation schemes on hanging nodes were applied. A cell matching method was used to handle periodic boundary condition. The difference schemes at nonuniform mesh for the universal Reynolds equation were derived. Ausas' cavitation algorithm was integrated into the AMR algorithm for the calculation of the bearing cavitation solution. The Richardson extrapolation method was employed as an a posteriori error estimation to tag the areas where they need to be refined. The AMR procedure was programmed in MATLAB. The details of the mesh management method, main algorithmic logics, nonuniform difference schemes, and error estimation were presented below.

Mesh and Data Storage
The example of a cell-based mesh is shown in Figure 2. The mesh is based on a quadtree data structure: each child cell owns a parent cell, and each parent cell possesses four child cells. The mesh structure is managed and stored in a matrix table, as shown in Table 1. With the matrix table, the operations on refining a leaf cell into four new leaf cells and coarsening four leaf cells back into their parent cell are easily implemented. The leaf cells mean the child cells without descendants and are the cells that actually participated in numerical calculation. Additionally, query operations can be applied within the matrix table to trace the ancestors and descendants for any cell according to need. More detailed quadtree data structure can be seen in [33]. The Mesh refining () and Mesh coarsening () algorithms are shown in Algorithms 1 and 2, respectively In the paper, a cell-based AMR algorithm was developed for bearing cavitation simulation. The aim is twofold: (1) The AMR method is introduced into cavitation problem. It provides another way to significantly improve the computational efficiency. This method exhibits potential advantages in analysis of bearings with complex geometric shapes. (2) Due to the complexity of the AMR method, the basic AMR algorithm was explained in the paper. The algorithm may provide guidance for people interested in AMR. The AMR algorithm was constructed based on a quadtree data structure with a Z-order filling curve managing cells. The hybrids of interpolation schemes on hanging nodes were applied. A cell matching method was used to handle periodic boundary condition. The difference schemes at nonuniform mesh for the universal Reynolds equation were derived. Ausas' cavitation algorithm was integrated into the AMR algorithm for the calculation of the bearing cavitation solution. The Richardson extrapolation method was employed as an a posteriori error estimation to tag the areas where they need to be refined. The AMR procedure was programmed in MATLAB. The details of the mesh management method, main algorithmic logics, nonuniform difference schemes, and error estimation were presented below.

Mesh and Data Storage
The example of a cell-based mesh is shown in Figure 2. The mesh is based on a quadtree data structure: each child cell owns a parent cell, and each parent cell possesses four child cells. The mesh structure is managed and stored in a matrix table, as shown in Table 1. With the matrix table, the operations on refining a leaf cell into four new leaf cells and coarsening four leaf cells back into their parent cell are easily implemented. The leaf cells mean the child cells without descendants and are the cells that actually participated in numerical calculation. Additionally, query operations can be applied within the matrix table to trace the ancestors and descendants for any cell according to need. More detailed quadtree data structure can be seen in [33]. The Mesh refining () and Mesh coarsening () algorithms are shown in Algorithms 1 and 2, respectively Algorithm 1 Mesh refining ()   for each cell  if a cell is a leaf cell and the refinement flag is 1  add four lines at the end of the matrix table; fill mesh relation information (cell number, type, parent, level, xy coordinate, and so on); interpolate initial flow variable using scatteredInterpolant() function (MATLAB function); end end add four lines at the end of the matrix table; fill mesh relation information (cell number, type, parent, level, xy coordinate, and so on); interpolate initial flow variable using scatteredInterpolant() function (MATLAB function); end end Algorithm 2 Mesh coarsening () for each cell if four leaf cells belong to a common parent cell and all coarse flags are 1 remove the lines of the four leaf cells; modify the parent cell into leaf cell; interpolate initial flow variable using the average mean of the four child cells; end end

Neighbor Finding
In a numerical solution of PDE, the calculation on partial derivatives needs to use the information of neighbor cells. It is necessary to find all the neighbor cell numbers for a specified cell. The essential way to find the neighbors is to perform a mirror query

Neighbor Finding
In a numerical solution of PDE, the calculation on partial derivatives needs to use the information of neighbor cells. It is necessary to find all the neighbor cell numbers for a specified cell. The essential way to find the neighbors is to perform a mirror query operation, as shown in Figure 3. The algorithm can be divided into two steps: find neighbor of same size and find neighbor of smaller size. A more detailed explanation is illustrated in [34,35]. The Neighbor finding () algorithm is shown in Algorithm 3.

Algorithm 3
Neighbor finding () find neighbor of same size using ADJ function [34] according to a given direction; store query path simultaneously; find neighbor smaller size according to mirror query path using REFLECT function [34]; Appl. Sci. 2021, 11, x FOR PEER REVIEW 5 of 15 operation, as shown in Figure 3. The algorithm can be divided into two steps: find neighbor of same size and find neighbor of smaller size. A more detailed explanation is illustrated in [34,35]. The Neighbor finding () algorithm is shown in Algorithm 3.

Algorithm 3
Neighbor finding () find neighbor of same size using ADJ function [34] according to a given direction; store query path simultaneously; find neighbor smaller size according to mirror query path using REFLECT function [34];

Neighbor Fixing
AMR procedure usually requires the mesh satisfying 2:1 condition between neighbors, meaning the number of neighboring cells to a specified cell is not allowed to exceed three. In other words, the level difference between adjacent cells should be limited to 0, 1, and −1. The 2:1 condition should be checked by the Neighbor finding () algorithm after mesh refinement and coarse operations. If the condition is not satisfied, the cells with more than three neighbor cells are divided into four cells by Mesh refining () algorithm operation.
The 2:1 condition in the diagonal direction is also mandatory in the present AMR procedure. If the condition is not implemented, the interpolation scheme on hanging nodes becomes much complex, which is not what we want. Figure 4 shows the mesh satisfying 2:1 condition after the Neighbor fixing () algorithm operation. The mesh was fixed from the mesh of Figure

Neighbor Fixing
AMR procedure usually requires the mesh satisfying 2:1 condition between neighbors, meaning the number of neighboring cells to a specified cell is not allowed to exceed three. In other words, the level difference between adjacent cells should be limited to 0, 1, and −1. The 2:1 condition should be checked by the Neighbor finding () algorithm after mesh refinement and coarse operations. If the condition is not satisfied, the cells with more than three neighbor cells are divided into four cells by Mesh refining () algorithm operation.
The 2:1 condition in the diagonal direction is also mandatory in the present AMR procedure. If the condition is not implemented, the interpolation scheme on hanging nodes becomes much complex, which is not what we want. Figure 4 shows the mesh satisfying 2:1 condition after the Neighbor fixing () algorithm operation. The mesh was fixed from the mesh of Figure 2. The 24 cells indicated by red lines are newly added. The Neighbor fixing () algorithm is shown in Algorithm 4.

Handingnode Interpolating
Hanging nodes are created at the interfaces between different level cells (see in Figure  2) [36]. Two cases need to be considered: (1) the calculation upon the coarse grids needs the interpolation of fine grids and (2) the calculation upon the fine grids needs the interpolation of coarse grids, as shown in Figure 5. In the case of (1), the value of the ghost node is simply calculated by the average value as where φ is the flow variable. In the case of (2), the value of the ghost node is calculated by quadratic interpolation as where a, b, c are the coefficients determined by φ2, φ3, φ4, and xi is the x coordinate of φ1.
The quadratic interpolation is implemented by MATLAB function [37]. If a boundary or fine grid is encountered, case (2) is considered a particular case. The values of φ2 and φ4 are, respectively, taken as the value of the boundary condition and the average value of the fine grid, as follows where φb is the value of the boundary condition. The present interpolation hybridizes linear and quadratic rules. Although all quadratic interpolations are feasible, more effort in programming is required. Other methods such as prolongation and restriction methods to handle interpolation can be seen in [38]. The Handingnode interpolating () algorithm is shown in Algorithm 5.

Handingnode Interpolating
Hanging nodes are created at the interfaces between different level cells (see in Figure 2) [36]. Two cases need to be considered: (1) the calculation upon the coarse grids needs the interpolation of fine grids and (2) the calculation upon the fine grids needs the interpolation of coarse grids, as shown in Figure 5. In the case of (1), the value of the ghost node is simply calculated by the average value as where ϕ is the flow variable. In the case of (2), the value of the ghost node is calculated by quadratic interpolation as where a, b, c are the coefficients determined by ϕ 2 , ϕ 3 , ϕ 4 , and x i is the x coordinate of ϕ 1 . The quadratic interpolation is implemented by MATLAB function [37]. If a boundary or fine grid is encountered, case (2) is considered a particular case. The values of ϕ 2 and ϕ 4 are, respectively, taken as the value of the boundary condition and the average value of the fine grid, as follows where ϕ b is the value of the boundary condition. The present interpolation hybridizes linear and quadratic rules. Although all quadratic interpolations are feasible, more effort in programming is required. Other methods such as prolongation and restriction methods to handle interpolation can be seen in [38]. The Handingnode interpolating () algorithm is shown in Algorithm 5.

Algorithm 5 Handingnode interpolating ()
if the neighbor cell is a boundary (Neighbor finding () algorithm returns zero) return the value of boundary condition; else if the number of neighbor cells is two (case (1)) return the average value using Equation (1); end if the number of neighbor cells is one and the level of neighbor cells is small than the level of the central cell (case (2)) return the quadratic interpolation value using Equation (2); end if the neighbor cell is the particular case (2) return the quadratic interpolation value using Equation (3); end end return the quadratic interpolation value using Equation (2); end if the neighbor cell is the particular case (2) return the quadratic interpolation value using Equation (3); end end

Periodic Matching
The west boundary and east boundary of the computational domain are paired into a periodic boundary condition. A periodic boundary condition is defined for the two boundaries where their values are linked in some defined way, such as interpolation or matching method. A matching method is applied to realize the periodic boundary condition. There are two steps to achieve the aim: (1) store the cell numbers of the west boundary and east boundary, respectively, and (2) conduct a pair matching operation according to the size of a cell.
The matching method is illustrated in Figure 6 and Table 2. In the first step, the cells of the west boundary and east boundary are detected by the Neighbor finding () algorithm according to the returned number of neighbor cells be zero or not. Owing to the Z-order filling curve (see below), the cell number is naturally in descending order along the y coordinate. Otherwise, a sort of operation needs to be performed by comparing y coordinate. Secondly, a pair number is added at each cell starting from the first two cells. The pair number is self-added when the size of the west boundary cell is equal to the size of the east boundary cell. For example, the cell of 28,95,97 matches the cell of 86 with the pair number of 2. Moreover, due to the number of neighbor cells of 86 exceeding three, dissatisfying the 2:1 condition, a mesh refinement operation should be conducted to the cell of 86 in the next step. The Periodic matching () algorithm is shown in Algorithm 6.

Periodic Matching
The west boundary and east boundary of the computational domain are paired into a periodic boundary condition. A periodic boundary condition is defined for the two boundaries where their values are linked in some defined way, such as interpolation or matching method. A matching method is applied to realize the periodic boundary condition. There are two steps to achieve the aim: (1) store the cell numbers of the west boundary and east boundary, respectively, and (2) conduct a pair matching operation according to the size of a cell.
The matching method is illustrated in Figure 6 and Table 2. In the first step, the cells of the west boundary and east boundary are detected by the Neighbor finding () algorithm according to the returned number of neighbor cells be zero or not. Owing to the Z-order filling curve (see below), the cell number is naturally in descending order along the y coordinate. Otherwise, a sort of operation needs to be performed by comparing y coordinate. Secondly, a pair number is added at each cell starting from the first two cells. The pair number is self-added when the size of the west boundary cell is equal to the size of the east boundary cell. For example, the cell of 28,95,97 matches the cell of 86 with the pair number of 2. Moreover, due to the number of neighbor cells of 86 exceeding three, dissatisfying the 2:1 condition, a mesh refinement operation should be conducted to the cell of 86 in the next step. The Periodic matching () algorithm is shown in Algorithm 6.

Z-Order Filling Curve
Z-order filling curve [39] is used to manage the solution order of the nonuniform mesh. A vector storing the ordered cell numbers is generated by the Z-ordering () algorithm, where the sequence of the cell numbers obeys the rule of the Z-order filling curve, as shown in Figure 7. The basic idea of the Z-ordering () algorithm is to travel the quadtree transversely under the sequence of NW→NE→SW→SE layer by layer. The Z-ordering () algorithm is shown in Algorithm 7.

Z-Order Filling Curve
Z-order filling curve [39] is used to manage the solution order of the nonuniform mesh. A vector storing the ordered cell numbers is generated by the Z-ordering () algorithm, where the sequence of the cell numbers obeys the rule of the Z-order filling curve, as shown in Figure 7. The basic idea of the Z-ordering () algorithm is to travel the quadtree transversely under the sequence of NW→NE→SW→SE layer by layer. The Zordering () algorithm is shown in Algorithm 7.

Difference Schemes on Nonuniform Mesh
The universal Reynolds equation [18] governing both liquid and cavitation regions is expressed as ∂ ∂x h 3 ∂p ∂x where h is the film thickness, p is the liquid pressure, µ is the dynamic viscosity, U is the sliding speed, and θ is the density ratio. As mentioned above, AMR allows different numerical resolutions to exist in the computational domain; therefore, the regions have different sizes. As a result, the difference scheme for nonuniform mesh is different from that for uniform mesh. Figure 8 shows the parameter definition for the nonuniform difference scheme.
where h is the film thickness, p is the liquid pressure, μ is the dynamic viscosity, U is the sliding speed, and θ is the density ratio. As mentioned above, AMR allows different numerical resolutions to exist in the computational domain; therefore, the regions have different sizes. As a result, the difference scheme for nonuniform mesh is different from that for uniform mesh. Figure 8 shows the parameter definition for the nonuniform difference scheme. The difference scheme of the pressure flow term is derived as 3  3  3  3  3  3  3  3  2  1  3   1  1  2  1 2  1  2  2  1  2   3  3  3  3  3  3  3  3  4  3  3   3 3 The difference scheme of the shear flow term is derived as Then, the explicit iterative schemes of p and θ can be written in the forms: where k is the number of iterations. The successive over-relaxation (SOR) method is applied to accelerate convergence. The SOR factors for p and θ are both taken as 1.2. The detailed solution strategy for Ausas' algorithm can be seen in [18,19].

Error Estimation
Discretization error is defined as the difference between the exact solution of the discrete equation and the exact solution of the PDE. It is the primary source of numerical The difference scheme of the pressure flow term is derived as The difference scheme of the shear flow term is derived as Then, the explicit iterative schemes of p and θ can be written in the forms: where k is the number of iterations. The successive over-relaxation (SOR) method is applied to accelerate convergence. The SOR factors for p and θ are both taken as 1.2. The detailed solution strategy for Ausas' algorithm can be seen in [18,19].

Error Estimation
Discretization error is defined as the difference between the exact solution of the discrete equation and the exact solution of the PDE. It is the primary source of numerical errors compared with round-off error and iterative error [40]. The Richardson extrapolation method is a recovery method and can be used to estimate discretization error [41].
The Richardson extrapolation method can obtain a more accurate solution than the available solution of the finest mesh. Assuming two numerical solutions have been calculated [42]: a fine grid with grid size h 1 and computed solution fc 1 ; and a coarse grid with grid size h 2 and computed solution fc 2 , the extrapolated solution is evaluated as where r is the grid refinement ratio r = h 2 /h 1 and p is the order of accuracy. In other words, the more accurate solution is obtained by fine grid solution f c1 and coarse grid solution f c2 with weight coefficients of r p /(r p − 1) and -1/(r p − 1). The discretization error for the fine grid solution f c1 could be estimated as follows In the AMR method, the present existing mesh is treated as the fine mesh. The coarse mesh is generated by coarsening the whole mesh: all leaf cells are replaced by their parent cells. The value of r is therefore ensured to be constant 2. The order of accuracy p is also treated as constant. Thus, the denominator of Equation (9) is unchanged. The error estimation is taken as Equation (10) indicates that the difference between the solutions obtained on the two meshes at each point is proportional to the local discretization error at that point [27].
When the error estimation operation is finished, the next step is to decide the value of the unacceptably large error. A practical way is to sort these errors from large to small, and to refine the top cells with large errors. The number of the cells need to be refined is determined according to need. The mesh refinement operation can be repeated until the solution is satisfactory.

Results and Discussion
The implementation steps of the AMR method are to (1) obtain a preliminary solution on initial coarse mesh, (2) obtain a more accurate solution by a refined mesh operation, and (3) repeat step (2) until a sufficiently accurate solution satisfying the need. The refinement strategy has a significant effect on the final mesh. It affects not only the accuracy of the solution but also the computational cost. An appropriate refinement strategy can realize a more efficient solution.
Although the use of variable gradient as a refinement indicator is more natural, the paper adopted a posteriori error estimation to refine the areas where the discretization error was large. As mentioned above, the Richardson extrapolation method is a simple way to find the area where it needs to be refined.
Pressure and density ratio are the main flow variables of a bearing film. If an individual parameter (pressure or density ratio) is used as error estimation, this is inappropriate in Ausas' cavitation algorithm since the pressure is always constant cavitation pressure in the cavitation region, and the density ratio is always constant 1 in the liquid region. The refinement on the whole bearing film is impossible. Significant discretization error always exists in the cavitation region or liquid region. Thus, pressure and density ratio are both taken as the variables to estimate error. Two cases were studied below to show the AMR results.
The first case is the cavitated journal bearing studied by Cupillard et al. [8,9] and Brewe [10]. The bearing is a finite type, with L/D = 4/3 (D = 100 mm), operating in the condition of ε = 0.60, n = 459 r/min, p c = 28 kPa (a), and µ = 0.0127 Pa·s. The solution of the three-time refined mesh is compared with the solution of the corresponding uniform mesh and the experimental data. The initial mesh is taken as a very coarse mesh with 16 × 16 cells. This mesh is generated by recursively refining the whole mesh four times. In the first refinement, 30 cells in the liquid region and 40 cells in the cavitation region were tagged to be refined. The tagged cells were the cells of the coarsen mesh with 8 × 8 cells in the implementation of the Richardson extrapolation method. After the refinement operation, the number of the cells refined was 240. It was lower than the expected number of grids, 280, as calculated by 70 × 4. This is because some of the tagged cells with large errors were common in the cavitation and liquid regions. In the subsequent refinements, 90 and 300 cells in the liquid region and 40 and 180 cells in the cavitation region were tagged to be refined for the second and third refinements, respectively. It is noted that the number of cells tagged was empirically determined according to the initial numerical tests. Although the refinement strategy is empirical, it still shows the law that the ratio of the tagged number to the total number reduces with the increase in total cell numbers. As shown in Table 3, the refinement ratios are 91%, 65%, and 41% in descending order. The refinement strategy indicates that the discretization error of the coarse mesh should be reduced sufficiently. This is because the coarse mesh has the largest discretization error, where the discretization error is proportional to the power p of the grid spacing h.  Figure 9 shows the pressure and density ratio for the final refined mesh. It can be seen that the refined meshes tagged by the pressure-based error estimation were basically distributed in the region with large pressure gradients, and the refined meshes tagged by density ratio-based error estimation were basically distributed in the region of cavitation interfaces where the density ratio gradients were also large. This behavior indicated that the discretization error was concentrated in the area with large variable gradients. Figure 10 compares the pressure distribution for the AMR results, uniform mesh results, and experimental results provided by Jakobsson and Floberg [12] and Olsson [13]. The pressure obtained by the AMR method is nearly consistent with that obtained by the uniform mesh (128 × 128 cells), and there is a slight difference in the experimental results. The number of the final refined mesh is 7984, while the number of the uniform mesh is 16,384. The AMR method provides almost the same accuracy results while the number of cells is reduced to around 50% of the uniform mesh. The computational efficiency is shown to be effectively improved.   The second case is the dimple-enhanced seal-like thrust bearing investigated by Qiu and Khonsari [22]. The bearing surface was textured with duplicated micro circular holes. The variation of the hole height results in a divergence and convergence clearance region. Hence, oil film rupture and reformation simultaneously occur around the hole. The bearing was operated in the condition of r0= 750 μm, hg = 10 μm, L = 3 mm, h0 = 4 μm, μ = 0.0035 Pa·s, U = 1.45 m/s, pc = 0.9 × 10 5 Pa (a), pa = 1.0 × 10 5 Pa (a). The refinement strategy for the bearing is listed in Table 3. The refinement ratios were 91%, 76%, and 48% for the three-time refinements. They are higher than those of case 1. A more computational cost was indicated to be paid. The solution of the three-times refined mesh is compared with the solution of the corresponding uniform mesh and the solution provided by Qiu and Khonsari [22], as shown in Figures 11 and 12. It can be seen that the AMR solution is close to the solution of the uniform mesh (128 × 128 cells) with a small difference compared with the solution provided by Qiu and Khonsari [22]. The difference is mainly reflected by the value of maximum pressure. Ausas' algorithm is sensitive to the number of grids. Insufficient number of grids often leads to high maximum pressure. The difference is acceptable under the current refinement strategy. The number of the final refined mesh is 9220, which is 56% of the uniform mesh. The AMR method provides almost the same accuracy solution while the computational cost is saved. The second case is the dimple-enhanced seal-like thrust bearing investigated by Qiu and Khonsari [22]. The bearing surface was textured with duplicated micro circular holes. The variation of the hole height results in a divergence and convergence clearance region. Hence, oil film rupture and reformation simultaneously occur around the hole. The bearing was operated in the condition of r 0 = 750 µm, h g = 10 µm, L = 3 mm, h 0 = 4 µm, µ = 0.0035 Pa·s, U = 1.45 m/s, p c = 0.9 × 10 5 Pa (a), p a = 1.0 × 10 5 Pa (a). The refinement strategy for the bearing is listed in Table 3. The refinement ratios were 91%, 76%, and 48% for the three-time refinements. They are higher than those of case 1. A more computational cost was indicated to be paid. The solution of the three-times refined mesh is compared with the solution of the corresponding uniform mesh and the solution provided by Qiu and Khonsari [22], as shown in Figures 11 and 12. It can be seen that the AMR solution is close to the solution of the uniform mesh (128 × 128 cells) with a small difference compared with the solution provided by Qiu and Khonsari [22]. The difference is mainly reflected by the value of maximum pressure. Ausas' algorithm is sensitive to the number of grids. Insufficient number of grids often leads to high maximum pressure. The difference is acceptable under the current refinement strategy. The number of the final refined mesh is 9220, which is 56% of the uniform mesh. The AMR method provides almost the same accuracy solution while the computational cost is saved.
(a) (b) Figure 11. Final refined mesh for thrust bearing: (a) pressure distribution; (b) density ratio distribution. Figure 12. Comparison of pressure distribution between different results for y = 1.5 mm.

Conclusions
A quadtree-based AMR procedure was developed for bearing cavitation simulation. The universal Reynolds equation was numerically solved on the nonuniform mesh by Ausas' cavitation algorithm. The pressure and density ratio were both taken as the indicators of the error estimation with the Richardson extrapolation method. The investigated cases showed that the AMR method provided nearly the same accuracy results compared with the uniform mesh, while the number of mesh was reduced to 50-60% of the number of the uniform mesh. The computational efficiency was effectively improved. The AMR method is suggested to be a potential tool for bearing cavitation simulation.

Conclusions
A quadtree-based AMR procedure was developed for bearing cavitation simulation. The universal Reynolds equation was numerically solved on the nonuniform mesh by Ausas' cavitation algorithm. The pressure and density ratio were both taken as the indicators of the error estimation with the Richardson extrapolation method. The investigated cases showed that the AMR method provided nearly the same accuracy results compared with the uniform mesh, while the number of mesh was reduced to 50-60% of the number of the uniform mesh. The computational efficiency was effectively improved. The AMR method is suggested to be a potential tool for bearing cavitation simulation.