Numerical calculation of oil film for ship stern bearing based on matrix method

Radial sliding bearings are widely used in ship shafting, its characteristics of lubricating oil film have important influence on the normal operation of the whole shaft system. In this work, the difference equations which is used to calculate the radial sliding bearing oil film features is transformed into matrix equations, the solving process be converted into solving matrix equation, combined with the powerful matrix calculation function of MATLAB, the solution process is simplified. It is not necessary to set the error precision and relaxation factor, so as to avoid the problem that the calculation result is not stable or even not convergent in the process of Successive Over Relaxation(SOR) method, and the calculation precision and stability are improved. The numerical results of matrix calculation method is compared with the result of SOR method, verified the correctness and feasibility of the matrix calculation method. Because the calculation is relatively stable, the matrix calculation method is more suitable for the calculation core of the relative computing software.


Introduction
Oil-lubricated radial sliding bearing is widely used in ship propulsion shafting, and its good lubrication and stability is related to the safety and normal operation of the whole ship. Therefore, the research on the lubrication performance of radial sliding bearing has been paid great attention to by scientific research institutions and scholars all over the world [1]. Nowadays, With the development of shipping industry and ships become high speed and large-scale, the requirements for bearings increased. At the same time, the development of modern computer technology also powerfully promoted the related researching works. In recent years, a large number of numerical calculation methods have been proposed for the lubrication characteristics of oil lubricated bearings, but summarize them, they are all in the process of solving the two-dimensional Reynolds equation [2].
In the study of fluid lubrication, the core is solving the Reynolds equation. But for the finite length bearings, the Reynolds equation is not analytic. Therefore, for a long time, researchers have obtained the approximate analytic solution through a large number of simplifying and assumptions, so as to obtain the bearing correlation characteristic parameters [3]. D 'agostino V et al adopted infinite long bearing hypothesis to obtain approximate analytical solutions for the Reynolds equation of finite length porous bearings [4]; Vignolo G takes the square of the aspect ratio: length over diameter(L/D)2, and the eccentricity ratio as perturbation and expansible parameters, proposes an analytical approximate solution of the Reynolds equation for isothermal finite length journal bearings by means of the regular perturbation method. When ( L/D)2 and the eccentricity ratio is small ,the method has high precision [5]. Sfyris D adopts zero pressure boundary and ignores the influence of eddy cutting velocity, and the approximate analytical solution of Reynolds equation is obtained by using the separation variable method and power series method [6]. Chasalevris A et al established the Reynolds equation for the three-oil leaf bearing and obtained its approximate solution by using the power series method [7]; Zhang YF et al. adopted the Reynolds boundary condition, used the separation variable method to solve the problem, and compared the results with the finite element simulation results [8]. Xiao Zhonghui et al. used a variational inequality to solve the oil film force and transformed the Jacobi matrix of the oil film into the linear equations, which reduced the calculation workload in a certain extent [9]. Lu Yanjun et al. used the finite element method to solve the oil film force and the Jacobi matrix by using the variational constraint principle [10]. H Di established a database for the oil film force of single wafers bearings, and use the interpolation operation to solve nonlinear oil film force [11]. Meng Zhiqiang based on the semi-Sommerfeld boundary conditions, using the separation of variables method to solve the finite-bearing nonlinear oil film force [12]. In general, the numerical solution method for Reynolds equation has been improved greatly in solving methods and solving precision, but the computational workload is still relatively large.
During the process of solving Reynolds equation, it is very important to select the appropriate calculation algorithm to increase the calculate precision and reduce the workload [13]. At present, there are many numerical methods for solving Reynolds equation. The most commonly used methods are finite difference method, finite element method and boundary element method. Finite difference method is relatively simple and it is widely used [14,15]. Usually the Reynolds equation is discretized by finite difference method and solved by Iterative solution method. In actual calculation, due to the nonlinear nature of the equation, the calculation results may not stable enough, and even not converge [14][15][16][17][18].
The present work based on the finite difference method, transformed The discretized difference equation into a matrix equation in a fixed format. Taking advantage of the powerful matrix calculation function of MATLAB to solve the Reynolds equation under Gumbel boundary conditions (Semi-Sommerfeld conditions), and the related parameters of the lubricating oil film of radial sliding bearing are obtained. The matrix calculation method is no need to set the relaxation factor and error precision, effectively avoid the disadvantage of the traditional SOR algorithm. It provides a relative stable way to solve the Reynolds equation.

Matrix solution based on the Gumbel conditions
According to the basic theory of bearing lubrication, under medium and low speed, medium and high load conditions, Reynolds equation of oil-lubricated radial sliding bearing [14] can be expressed as: 3 3 1 5 Define the top of the cross section of the bearing as 0 degree, The solution area of bearing is divided into m and n equidistant grids along the bearing circumference and axial direction respectively, as is shown in figure 1. The position of 0°and 360°are overlap, so the circumferential is divided into m nodes, axial divided into n+1 node. The bearing is simplified into m (n+1) nodes. The oil film pressure and thickness of any node Ai, j are expressed by Pi, j and hi,j . Define the environmental pressure P=Pa, The bearing length is L. The terminal pressure of bearing equals environmental pressure, That is : when j = 1& n+1, z = 0& L, Pi,1=Pi,n+1= Pa; Therefore, there just remaining m(n-1) nodes that the pressure values need to be calculated. For any node Ai,j, according to the difference principle, the first-order and second-order partial derivatives can be expressed by the four pressure values of the nodes around it, as is shown in figure 2. According to the relationship, Formula 1 can be converted to into the form as follow. 2  2   3  3  2  3  , , , Regardless of the Angle between shaft and the bearing, the journal center line is parallel to the center line of the bearing, the oil film thickness can be calculated easily [14]. And the oil film thickness keep the same in the axial direction, that is: , for each node, a difference equation can be built. Totally m(n-1) independent equations can be listed for the whole simplified model. Write equation 2 in a another way.
In which: In order to facilitate the simplification of matrix equation and the establishment of matrix model, each node is numbered, as is shown in figure 3-a. The grid is processed:

Calculation example
For the sake of easy understanding, take the grid parameter m=4, n=3 for example to establish the matrix equation. Assuming that the pressure of the eight nodes on the two ends of the bearing, z=0 and z=L, is equal to the environmental pressure Pa; The corresponding difference equations are built for the unknown 8 nodes as follow.
Each units In the matrix equation can be written:

SOR algorithm based on the Gumbel condition
In 1970s, D. M. Young proposed Successive Over Relaxation method, called SOR method for short, is one of the effective methods to solve large sparse matrix equations [14,15]. The key to use the over relaxation iteration method is to choose the appropriate relaxation factor. If the relaxation factor is selected properly, it will greatly shorten the computation time. The deformation of equation 6 can be written as: In which, ω is the relaxation factor. The key to solve Reynolds equation with SOR algorithm is the selection of error accuracy and relaxation factor. Take  =1.0  10-5, ω=1.9, calculate the stern bearing in the third section above with the SOR algorithm. The result of distribution of oil film pressure in the bearing is shown in figure 5.

Comparison of the two solutions
In order to verify the reliability of the method, two calculation methods are used to calculate the actual bearing in section 3 based on Gumbel boundary conditions, and the results are compared. As a common calculation method, the calculation results of SOR algorithm are considered to be correct and reliable in this paper. The environmental pressure is defined : Pa=0.1MPa. Figure. 6(a) and (b) and Figure.   It can be seen from figure 6 and figure 7, the calculation results of using the matrix method and SOR algorithm are basically identical. The maximum pressure value in both calculation are occurred in 156°circumferential position, and the maximum value is 5.651 MPa. The oil film pressure drops to 0 MPa at 180°circumferential position, which satisfies the Gumbel condition. Both calculations do not consider the angle between shaft and bearing axis. Therefore, the pressure of the oil film is symmetrical distribution along the axial direction of the bearing. The highest pressure is located in the center of the bearing, and the oil pressure on both sides of the bearing is equal to the pressure of the environment. Although the maximum pressure value is consistent, there are small differences in the oil pressure distribution because of the difference of two calculation methods, relatively, the difference little and it can be neglected. The comparison results show that the calculation accuracy of matrix calculation method is completely satisfied with the actual calculation, and the calculation results are accurate and reliable.
When the SOR algorithm is used in calculation, the selection of the relaxation factor and the error accuracy value will have a great influence on the calculation precision and the calculation efficiency. In table 1, the influence of different relaxation factors and error accuracy values on the calculation results and calculation efficiency is presented in calculating the example above by the SOR method. In table 1, F stands for the oil film force, x_z location is the equivalent acting position on the z axis of the force in the x direction, y_z location is the equivalent acting position on the z axis of the force in the y direction, T means the calculation time of the process. In theory, the equivalent position of the oil film counterforce on the z axis should be in the middle of the bearing length, where L=530mm. Data in table 1 shows that the smaller the error accuracy, the closer the calculation result is to the theoretical value, but the number of iterations increases and the calculation time increases. With the increase of the relaxation factor, the calculation time decreases and then increases, and the calculation results are stable. When the relaxation factor is equal to 2.00, the result is divergent. The calculation results show that when relaxation factor equals 1.9,the calculation time is shorter and the calculation efficiency is relative high. In fact, for different bearing under different working conditions, the optimal relaxation factor is different.
For the problem, the matrix solutions based on the Gumbel boundary condition show certain advantages. Without the process of setting the optimal relaxation factor and error precision, the whole calculation time is only 0.09s. In the case of relative complex computing processes ,in order to satisfy the same precision, the SOR method usually requires to calculate hundreds or even thousands of times, when the efficiency advantage of matrix solution will be more obvious.
In conclusion, due to the advantages of matrix solution, effectively prevent the computation of SOR may appear unstable problem. In terms of computational efficiency and stability calculation is better than the SOR method, the matrix method is more suitable as the computing core for a general calculation software.

Conclusion
The calculation equation of the lubricating oil film is directly transformed into a matrix equation, and the solving process is simplified by using the powerful matrix computing ability of MATLAB. Calculate the same bearing model by using the matrix method and SOR method respectively based on Gumbel conditions. The results show that the result of two methods are consistent, verified the accuracy and feasibility of matrix calculation method. By using the matrix method ,there is no necessary to set the relaxation factor and the error precision make the calculation process simpler and the calculation result is more stable. The matrix method effectively avoids the problems that may arise when using traditional SOR computing methods. It has more advantages in the development of computer software for bearing lubrication film calculation.