SUMMARYUtilizing the mathematical programming technique of contact problems and Fast Multipole Boundary Element Method (FM-BEM) theory, Taylor Series Multipole BEM (TSMBEM) for solving 3D multibodies elastic frictional contact problems is presented based on the 'node-to-surface' contact model. The numerical method is able to analyze the large-scale contact problems using non-conforming or conforming boundary meshes. For multi-bodies contact problems, the computational costs of contact identification can be effectively reduced by the contact table. The mathematical programming method for TSMBEM is presented based on the interpolation approach of contact elements. The linear system arising from multipole BEM is solved by using Generalized Minimal RESidual algorithm GMRES(m) in conjunction with near-field preconditioning technique. The numerical experiments of flat punch and HC roll system are carried out by the newly developed multipole BEM. The numerical results are compared with the Hertz theoretical solutions, which illustrate the effectiveness and the validity of the numerical method for 3D multi-bodies elastic frictional contact analysis.