Unlike the Navier boundary condition, this paper deals with the case when the slip of a fluid along the wall may occur only when the shear stress attains certain bound which is given a priori and does not depend on the solution itself. The mathematical model of the velocity–pressure formulation with this type of threshold slip boundary condition is given by the so-called variational inequality of the second kind. For its discretization, we use P1-bubble/P1 mixed finite elements. The resulting algebraic problem leads to the minimization of a nondifferentiable energy function subject to linear equality constraints representing the discrete impermeability and incompressibility condition. To release the former one and to regularize the nonsmooth term characterizing the stick–slip behavior of the algebraic formulation, two additional vectors of Lagrange multipliers are introduced. Further, the velocity vector is eliminated, and the resulting minimization problem for a quadratic function depending on the dual variables (the discrete pressure and the normal and shear stress) is solved by the interior point type method which is briefly described. To justify the threshold model and to illustrate the efficiency of the proposed approach, three physically realistic problems are solved and the results are compared with the ones solving the Stokes problem with the Navier boundary condition.