We present numerical approximations of the 3D steady state Navier-Stokes equations in velocity-pressure formulation using trivariate splines of arbitrary degree d d and arbitrary smoothness r r with r > d r>d . Using functional arguments, we derive the discrete Navier-Stokes equations in terms of B B -coefficients of trivariate splines over a tetrahedral partition of any given polygonal domain. Smoothness conditions, boundary conditions and the divergence-free condition are enforced through Lagrange multipliers. The pressure is computed by solving a Poisson equation with Neumann boundary conditions. We have implemented this approach in MATLAB and present numerical evidence of the convergence rate as well as experiments on the lid driven cavity flow problem.