In this paper, explicit boundary-domain integral equations for evaluating velocity gradients are derived from the basic velocity integral equations. A free term is produced in the new strongly singular integral equation, which is not included in recent formulations using the complex variable differentiation method (CVDM) to compute velocity gradients (Int. J. Numer Meth. Fluids 2004; 45:463-484; Int. J. Numer Meth. Fluids 2005; 47:19-43). The strongly singular domain integrals involved in the new integral equations are accurately evaluated using the radial integration method (RIM). Considerable computational time for evaluating integrals of velocity gradients can be saved by using present formulation than using CVDM. The formulation derived in this paper together with those presented in reference (Int. J. Numer. Meth. Fluids 2004; 45:463-484) for 2D and in (Int. J. Numer Meth. Fluids 2005; 47:19-43) for 3D problems constitutes a complete boundary-domain integral equation system for solving full Navier-Stokes equations using primitive variables. Three numerical examples for steady incompressible viscous flow are given to validate the derived formulations. Copyright 2007 John Wiley & Sons, Ltd.