The thin, flexible, flat belt system is frequently employed for power transmission, conveyor systems, and various manufacturing process equipments, such as printers, elevators, and automated teller machines. In general, some critical problems can happen during operation, since a thin flexible belt on rollers can easily skew. This is especially serious for high-speed belt operations. In this skewed situation, due to possible relative roller misalignment, roller deflection, and/or belt non-uniformity itself, a thin flexible belt can move along the perpendicular direction to the travelling direction. In this investigation, an efficient analysis method for simulating skewing phenomena of thin flexible belt systems is presented using the complete integration of large deformable finite element methods and multibody dynamics. A 4-node shell element is used to describe a thin flexible belt and two rigid cylindrical bodies describe the roller. Each cylindrical body is fixed with a pin joint and the positions of the rollers are adjusted to impose a tension to the belt. For the interactions between the belt and the rollers, an efficient contact search algorithm based on a compliant contact force model is presented in this study. The contact search algorithm is composed of a pre-search and a detailed search to identify points of contact between the belt and the rollers. In the pre-search, a simple algorithm is used to identify lines on the belt that are close enough to the rollers that there might potentially be contact between them. The detailed search is performed by making efficient two-dimensional (2D) comparisons of the cross-section of the roller with a 2D projection of the belt lines identified in the pre-search. Afterwards, the contact force is generated by a compliant force model that uses a stick-slip frictional algorithm. The simulation model is validated against experimental measurement systems for thin flexible belt and rollers. It shows a good agreement with the numerical results. Furthermore CPU time with commercial software ANSYS Workbench is compared to check the performance of the algorithm. The proposed method is more than 10 times faster compared with ANSYS Workbench.