Summary: In this work we propose a new stabilized approach for solving the incompressible Navier-Stokes equations on fixed overlapping grids. This new approach is based on the partition of unity finite element method, which defines the solution fields as weighted sums of local fields, supported by the different grids. Here, the discrete weak formulation of the problem is re-set in cG(1)cG(1) stabilized form, which has the dual benefit of lowering grid resolution requirements for convection dominated flows and allowing for the use of velocity and pressure discretizations which do not satisfy the inf-sup condition. Additionally, we provide an outline of our implementation within an existing distributed parallel application and identify four key options to improve the code efficiency namely: the use of cache to store mapped quadrature points and basis function gradients, the intersection volume splitting algorithm, the use of lower order quadrature schemes, and tuning the partition weight associated with the interface elements. The new method is shown to have comparable accuracy to the single mesh boundary-fitted version of the same stabilized solver based on three transient flow tests including both 2D and 3D settings, as well as low and moderate Reynolds number flow conditions. Moreover, we demonstrate how the four implementation options have a synergistic effect lowering the residual assembly time by an order of magnitude compared to a naive implementation, and showing good load balancing properties.