## 2 Computational Model
### Numerical Solution  
When dealing with EHL, it is crucial to adequately load the contact surfaces, as the distribution of hydrodynamic pressure, shear stress, contact pressure and friction of the asperity contact have a significant influence on the elastic deformation of the contact bodies and vice versa. Due to this, the force interaction between the contact bodies is simulated by the modal force of Equation (5), where the forces are applied to the boundary nodes on the contact surfaces. Depending on the lubrication regime, two solutions to hydrodynamic lubrication, the universal Reynolds Equation (7) and Average Reynolds Equation (11), were presented in the previous section. The Average Reynolds equation is simply the universal Reynolds equation augmented by the flow factors. Since the flow factors compare the flow in a rough bearing with that of a smooth bearing, their effect disappears with the thick film thickness. This effect can be observed in the dependence of the flow factors on the oil film thickness available in Ref. [32]. Therefore, the Average Reynolds equation is more general and can be used in both partial lubrication regime and hydrodynamic lubrication regime. In this paper, only the Average Reynolds equation is used. In practical terms, flow factor data may not be available for a contact pair of interest. In this case, it is more suitable to begin with the universal Reynolds equation. A simpler model imitates the crawl-walk-run approach and is easier to debug. Then the model should be refined using the Average Reynolds equation. The Reynolds equation (in both forms) can be further simplified for a specific application of interest. First, the lubricant density can be considered constant with a certain level of simplification [37; 38]. In this way, the hydrodynamic pressure does not depend on the lubricant density value. In addition, further simplification can be achieved if some surfaces move in one direction, and the surface velocity in other directions are neglected. If one of the bearing dimensions is significantly smaller than the other, the pressure gradient (the Poiseuille term) in the large dimension direction is negligible. This is the essence of the solution of an infinitely short and infinitely long bearing. In this case, it is possible to derive an analytic expression of the hydrodynamic pressure for these simplified bearings. If a finite length bearing is investigated, the Reynolds equation is solved numerically. It is classified as a nonlinear partial differential equation. These types of equations may be discretized by either the finite difference method, the finite volume method or the finite element method. In the case of this manuscript, the finite difference method is used--each derivative is replaced with an approximate difference formula at each nodal point, and the new set of equations is generated. There are three approaches to the solution of this new set of equations: the relaxation (iterative) methods, the direct matrix methods and the rapid methods (Fourier methods, for example) [39]. In the case of this manuscript, the Gauss-Seidel iteration method with the overrelaxation parameter \(\omega\) is used. The solution of the discrete node \(i,j\) is as follows:  
\[P_{i,j}^{(k+1)}=P_{i,j}^{(k)}-\omega\frac{r_{i,j}^{(k)}}{J_{i,j}}, \tag{19}\]  
where \(P^{(k)}\) indicates the dimensionless pressure of the previous iteration \((k)\), \(r\) reflects the residuum, and \(J\) represents the Jacobian. Moreover, the Chebyshev acceleration can be used for the overrelaxation parameter estimation [39].