With the increasing research in the field of contact mechanics, different types of contact models have been investigated by many researchers by employing various complex material models. To ascertain the orthotropy effect and modeling parameters on a receding contact model, the double frictional receding contact problem for an orthotropic bilayer loaded by a cylindrical punch is taken into account in this study. Assuming plane strain sliding conditions, the governing equations are found analytically using Fourier integral transformation technique. Then, the resulting singular integral equations are solved numerically using an iterative method. The weight function describing the asymptotic behavior of the stresses are investigated in detail and powers of the stress singularities are provided. To control the trustworthiness and correctness of the analytical formulation and to compare the resulting stress distributions and contact boundaries, a numerically efficient finite element method was employed using augmented Lagrange contact algorithm. The aim of this paper is to investigate the orthotropy effect, modeling parameters and coefficients of friction on the surface and interface stresses, surface and interface contact boundaries, powers of stress singularities, weight function and to provide highly parametric benchmark results for tribological community in designing wear resistant systems.