The Navier–Stokes equations play a key role in the modeling of blood flows in the vascular sys- tem. The cost for solving the 3D linear system obtained by Finite Element (FE) discretization of the equations, using tetrahedral unstructured meshes and time advancing finite difference schemes, is very high; to lower the time to solution and to address complex problems, using a parallel framework and developing specific preconditioners become necessary. Important factors to measure parallel performances of a preconditioner are the independence of the number of iterations with respect to the number of parallel tasks (scalability of the precondi- tioner), on the mesh size (optimality), and on the physical parameters (robustness), as well as its properties of strong and weak scalability. We propose a model to explain the effect of nonscalable preconditioners on a parallel iterative solver. Then we propose approximate versions of state of the art preconditioners for the Navier– Stokes equations, namely the Pressure Convection–Diffusion (PCD), Yosida, and SIMPLE preconditioners. We exploit factorizations of the linearized system where inverses are handled using specific core preconditioners such as, e.g., the algebraic additive Schwarz or algebraic multigrid preconditioners. We present new results for the Relaxed Dimensional Factorization (RDF) preconditioner which allows for an automatic parameter tuning. Weak and strong scalability results illustrate the efficiency of our approach using both classical benchmark problems and test cases relevant to hemodynamics simulations, using up to 8192 cores. Then, we extend our preconditioners for the Navier–Stokes equations to Fluid–Structure Interaction problems; we devise preconditioners that rely on physics–specific ad–hoc precon- ditioners for the fluid, the structure and the geometry subproblems. We compare the evolution of the number of iterations to solve the full system with classical methods on a geometry of physiological interest. We also investigate the strong scalability of our FSI solver. Finally, we describe how to develop a C++ flexible Navier–Stokes solver using a policy–based design. In particular we discuss different techniques to access a host class member variables from its compound policies and discuss their main advantages.