A spacecraft is a typical rigid–flexible–thermal coupled multibody system, and the study of such rigid–flexible–thermal coupled systems has important engineering significance. The dissipation effect of material damping has a significant impact on the response of multibody system dynamics. Owing to the increasing multitude of computational dimensions, computational efficiency has remained a significant bottleneck hindering their practical applications in engineering. However, due to the fact that the stiffness matrix is a highly nonlinear function of generalized coordinates, traditional methods of modal truncation are difficult to apply directly. In this study, the absolute nodal coordinate formulation (ANCF) is used to uniformly describe the modeling of rigid–flexible–thermal coupled multibody systems with large-scale motion and deformation. The constant tangent stiffness matrix and damping matrix can be obtained by locally linearizing the dynamic equation and heat transfer equations, which are based on the Taylor expansion. The dynamic and heat transfer equations obtained by reducing the order of complex modes are transformed into a unified first-order equation, which is solved simultaneously. The orthogonal complement matrix of the constraint equation is proposed to eliminate the nonlinear constraints. A strategy based on energy preservation was proposed to update the reduced-order basis vectors, which improved the calculation accuracy and efficiency. Finally, a systematic method for rigid–flexible–thermal coupled viscoelastic multibody systems via modal truncation with complex global modes is developed.