## 3 Numerical Solution and Model Validity  
### Numerical Solution of the System Model  
By combining the shaft model and quasi-static model of ACBB, the generic shaft-bearing system model is obtained. The unknowns in the system model include the ball orbital speed \(\omega_{c}\), spinning speed \(\omega_{b}\), pitch angle \(\beta\), axial and radial distances \((X_{1},X_{2})\) to imply the final position of the ball center, friction coefficients \(\mu\) at ball-raceway contacts, inner ring displacements \(\big{(}\delta_{x},\delta_{y},\delta_{z},\theta_{x},\theta_{y}\big{)}\) for ACBB, and the node displacements for the shaft. The unknowns are \(6Z+5\) for each ball bearing and \(5(N+1)\) for the shaft. Considering that the inner ring displacements of ball bearings are the same as the displacements of shaft nodes where ball bearings are mounted, the total unknowns of the shaft-bearing system model are \(k_{b}(6Z)+5(N+1)\), and \(k_{b}\) is the number of ACBBs supporting the shaft. Obviously, the shaft-bearing system model is a huge system of equations provided that the shaft is supported by several ACBBs and the number of ball elements for each bearing and of the beam elements for the shaft is large. It is computationally infeasible to solve all the equations of the shaft-bearing system model together. Hence, a numerical scheme is proposed to solve the system model. The numerical solution process is presented in Figure 6. The shaft-bearing system properties, including material properties, geometrical parameters, operation conditions, load conditions, axial preload, and/or angular misalignment of ball bearing, are taken as the input parameters. The initial guess for the node displacements of the shaft, the spinning speed \(\omega_{b}\), orbital speed \(\omega_{c}\), and pitch angle \(\beta\) of each ball element are also necessary to start the solution. Then, the distances parameters \((X_{1},X_{2})\) and friction coefficients \(\mu\) of all ball elements in the \(k\)th ball bearing are solved together by the global optimization method. The new spinning speed \(\omega_{b}^{*}\), orbital speed \(\omega_{c}^{*}\), and pitch angle \(\beta^{*}\) for ball elements are calculated based on Equations (1)-(3). The distances parameters \((X_{1},X_{2})\), friction coefficients \(\mu\), speed parameters (\(\omega_{b}\),\(\omega_{c}\)), and pitch angle \(\beta\) are solved alternately until the corresponding converge criterion is reached. Once the abovecalculations for all the ball bearings are finished, the equilibrium equations of the shaft are checked. Since the loads applied on the shaft nodes provided by the ball bearings are a nonlinear function of the node displacements of the shaft, the equilibrium equation of the shaft, Equation (25), is nonlinear as a function of the shaft nodes' displacements. Based on the first order Taylor expansion of nonlinear equilibrium equations of the shaft, the node displacements on the shaft are updated. The whole calculation process will be continued until the equilibrium equations of the shaft are satisfied. Then, the node displacements of the shaft, ball-raceway contact forces, spinning speed, orbital speed, and pitch angle for each ball element can be obtained. It should be noted that, in the innermost loop of the flowchart, the Newton-Raphson algorithm is not used to solve the equilibrium Equation (13) of all ball elements simultaneously to obtain the distance parameters \(X_{1},X_{2}\) and friction coefficients \(\mu\) for the \(k\)th ball bearing, although this method was commonly used in the previously published works. The accuracy of the Newton-Raphson algorithm typically relies on the trial-and-error of initial estimates, and the numerical scheme will be very time-consuming and not be able to converge if the initial chosen solution is far away from the exact solution. In order to overcome this deficiency, the global optimization method is used to calculate the unknowns \(X_{1},X_{2}\) and \(\mu\) of all ball elements simultaneously for the \(k\)th ball bearing.