In this contribution, we introduce a varying-order (VO) NURBS discretization method to enhance the performance of the isogeometric analysis (IGA) technique for solving three-dimensional (3D) large deformation frictional contact problems involving two deformable bodies. Building on the promising results obtained from the previous work on the 2D isogeometric contact analysis (Agrawal and Gautam, 2020), this work extends the method's capability for tri-variate NURBS-based discretization. The proposed method allows for independent, user-defined application of higher-order NURBS functions to discretize the contact surface while employing the minimum order NURBS for the remaining volume of the elastic solid. This flexible strategy enables the possibility to refine a NURBS-constructed solid at a fixed mesh with the controllable order elevation-based approach while preserving the original volume parametrization. The advantages of the method are twofold. First, employing higher-order NURBS for contact integral evaluations considerably enhances the accuracy of the contact responses at a fixed mesh, fully exploiting the advantage of higher-order NURBS specifically for contact computations. Second, the minimum order NURBS for the computations in the remaining bulk volume substantially reduces the computational cost inherently associated with the standard uniform order NURBS-based isogeometric contact analyses. The capabilities of the proposed method are demonstrated using various contact problems between elastic solids with or without considering friction. The results with the standard uniform order of tri-variate NURBS-based discretizations are also included to provide a comprehensive comparative assessment. We show that to attain results of similar accuracy, the varying-order NURBS discretization uses a much coarser mesh resolution than the standard uniform-order NURBS-based discretization, hence leading to a major gain in computational efficiency for isogeometric contact analysis. The convergence study demonstrates the consistent performance of the method for efficient IGA of 3D frictional contact problems. Furthermore, the simplicity of the method facilitates its direct integration into the existing 3D NURBS-based IGA framework with only a few minor modifications.