Contact problem simulation creates a significant challenge because of the nonlinearities and non-smooth characteristics of contact interaction. For contact modeling, isogeometric analysis (IGA) overcomes the limitations of finite element analysis (FEA) through the utilization of non-uniform rational B-spline (NURBS) functions. In contact simulation, IGA has demonstrated enhanced capabilities, outperforming FEA in robustness and efficiency. As contact is an interface phenomenon, efficiency of contact analysis depends on accurately representing contact surfaces. Within the IGA framework, the robustness and efficiency of a contact formulation can be improved through contact surface refinement using higher order NURBS functions. In this direction, several studies have been conducted, and a recent development is the varying-order (VO) NURBS discretization strategy. It has shown improved accuracy and reduced computer time compared to conventional NURBS discretization strategies. However, its application to frictional large deformation contact problems in three-dimensional (3D) settings has yet to be explored, which sets the goal of the current work. Herein, implementation for VO NURBS discretization is verified, showing the advantages over conventional NURBS discretization strategies in terms of accuracy and computer time.