In order to avoid the degeneration of high-speed spindle's machining accuracy in actual machining caused by the uneven distribution of temperature field at the design stage, a three-dimensional (3D) finite element analysis (FEA) model, which considered the combined influence of thermal contact resistance (TCR) and bearing stiffness on the accuracy of simulation results, was proposed to conduct transient thermal-structure interactive analysis of motorized spindles. And the method to calculate the boundary conditions used in the FEM model were discussed in detail, such as the heat loads, convective heat transfer, TCR and bearing stiffness. Based on the quasi-static mechanics analysis of rolling bearing, the transfer relationships among multiple variables in the equilibrium equation set of bearings were analyzed. Newton–Raphson algorithm, which regarded the contact angle as the iteration variable, was proposed to calculate the heat power and stiffness of bearings and improve the convergence of the algorithm, and the termination criteria of contact angle's searching trial calculation was proposed to improve the accuracy of the algorithm. The Weierstrass–Mandelbrot (W–M) function in fractal geometry was used to characterize the rough surface morphology of bearing rings. The fractal parameters were identified by the power spectrum method, and a contact deformation model of asperities was developed to calculate the contact parameters used in the TCR modeling. Then, the predictive model for TCR, which considered the combined effect of the morphology of bearing rings and the contact deformation of asperities, was proposed. Thermal equilibrium experiments were conducted to demonstrate the validity of the model. The results showed that the FEA model can accurately simulate the temperature field and thermal deformation of the spindle system and that the FEA model was much more accurate than the traditional thermal model of the high-speed spindle system which ignored TCR and bearing stiffness.