In this paper, three implementations of B-differentiable Newton method (BDNM) for solving 3D elastoplastic frictional contact problems are presented. In the first two methods, a nested strategy is employed. The contact equations, which are solely associated with the contact conditions in the normal and tangential directions due to the introduction of the contact flexibility matrix, are expressed as B-differentiable equations. They are solved in the external iteration, or regarded as a sub-step nested within the solution of nonlinear equilibrium equations, leading to the decoupling of the calculations for nonlinear equilibrium equations and contact forces. In the third method, both the nonlinear equilibrium equations and contact conditions are expressed as B-differentiable equations simultaneously. Both nodal displacement and contact force are treated as independent variables, solved using BDNM to ensure the global convergence, with a variant of a Jacobian matrix associated with unknowns derived. In the numerical examples, the accuracy of the proposed methods is demonstrated by comparing with the results obtained from ANSYS. The convergence and computational efficiency are investigated through a series of models with varying degrees of freedom (DOFs) and contact node pairs (CNPs). The robustness of the proposed methods is presented by applying them to the analysis of an arch dam with seven transverse joints.