This work presents a rigorous three-dimensional finite element procedure to analyze belt transmission systems. The frictional contact behavior between the belt and the pulley, which accounts for the power loss of the system and the wear of the belt, is investigated in detail. In addition to adopting the transformation matrix to satisfy the geometric conditions on the contact surfaces, the proposed procedure also uses the modified elements with incremental Wilson displacement modes to improve the accuracy due to bending at the end zones of the contact area for the belt. To demonstrate the accuracy and feasibility of the proposed procedure, the analyses for flat and V belt drives are carried out. Excellent correlations between the calculated results and referenced theoretical/experimental solutions are found. The influences of friction coefficients on the deformation, normal and tangential contact forces on the contact surfaces are studied as well. Those will be helpful for the estimation of wear properties and operation efficiency for belt transmission systems.