Axially symmetric contacts can be simulated rapidly with the Method of Dimensionality Reduction (MDR). Within its range of application, the method is highly accurate and several orders of magnitude faster than conventional finite element, or even boundary element approaches. The low computational complexity of the method is due to the low number of degrees of freedom, and it is due to the fact that these degrees of freedom are decoupled. The method was first proposed by V. L. Popov and S. Psakhie in 2007.

From a numerical perspective, the integral transformations which have to be performed as part of this method are a challenge. In the following paper, a simple method for the implementation of these transformations is proposed:

Fast numerical implementation of the MDR transformations
J. Benad
Facta Universitatis, Series: Mechanical Engineering, vol. 16, no. 2, pp. 127-138, 2018

The implementation technique has also been included in a book about the MDR by V. L. Popov, M. Heß and E. Willert.

It can be found here on ResearchGate.