A second-order accurate kernel-free boundary integral method is presented for Stokes and Navier boundary value problems on three-dimensional irregular domains. It solves equations in the framework of boundary integral equations, whose corresponding discrete forms are well-conditioned and solved by the GMRES method. A notable feature of this approach is that the boundary or volume integrals encountered in BIEs are indirectly evaluated by a Cartesian grid-based method, which includes discretizing corresponding simple interface problems with a MAC scheme, correcting discrete linear systems to reduce large local truncation errors near the interface, solving the modified system by a CG method together with an FFT-based Poisson solver. No extra work or special quadratures are required to deal with singular or hyper-singular boundary integrals and the dependence on the analytical expressions of Green's functions for the integral kernels is completely eliminated. Numerical results are given to demonstrate the efficiency and accuracy of the Cartesian grid-based method.