The particularities of the aircraft parts riveting process simulation necessitate the solution of a large amount of contact problems. We propose a primal-dual interior point method-based solver for solving such problems efficiently. The proposed method features a worst case polynomial complexity bound $O(\sqrt{n}\ln{\frac{1}{\epsilon}})$ on the number of iterations, where $n$ is the dimension of the problem and $\epsilon$ is a threshold related to desired accuracy. In practice, the convergence is often faster than this worst case bound, which makes the method applicable to large-scale problems. The computational challenge is solving the system of linear equations because the associated matrix is ill-conditioned. To that end, we introduce a preconditioner and a strategy for determining effective initial guesses based on the physics of the problem. We compare numerical results to ones obtained using the Goldfarb-Idnani algorithm. The results demonstrate the efficiency of the proposed method.