The main goal of this work is to devise robust iterative strategies to partition the solution of the Navier-Stokes equations in a three-dimensional (3D) domain, into non overlapping 3D subdomains, which communicate through the exchange of averaged/integrated quantities across the interfaces. The novel aspect of the present approach is that at coupling boundaries the conservation of flow rates and of the associated dual variables is implicitly imposed, entailing a weak physical coupling. For the solution of the non-linear interface problem two strategies are compared: relaxed fixed point and Newton iterations. The algorithm is tested in several configurations for problems ranging from academic ones to some related to the computational hemodynamics field, which involve more than two components at each coupling interface. In some cases it is shown that relaxed fixed point methods are not convergent, whereas the Newton method leads in all tested cases to convergent schemes. One of the appealing aspects of the strategy proposed here is the flexibility in the setting of boundary conditions at the interfaces, where no hierarchy is established a priori (unlike Gauss-Seidel methods). The usefulness of this methodology is also discussed in the context of dimensionally-heterogeneous coupling and preconditioning for domain decomposition methods.