We discuss aspects of implementation and performance of parallel iterative solution techniques applied to low Reynolds number flows around fixed and moving rigid bodies. The incompressible Navier-Stokes equations are discretised with Taylor-Hood finite elements in combination with a semi-implicit pressure-correction method. The resulting sequence of convection-diffusion and Poisson equations are solved with preconditioned Krylov subspace methods. To achieve overall scalability we consider new auxiliary algorithms for mesh handling and assembly of the system matrices. We compute the flow around a translating plate and a rotating insect wing to establish the scaling properties of the developed solver. The largest meshes have up to 132 × 10^6 hexahedral finite elements leading to around 3.3 × 10^9 unknowns. For the scalability runs the maximum core count is around 65.5 × 10^3. We find that almost perfect scaling can be achieved with a suitable Krylov subspace iterative method, like conjugate gradients or GMRES, and a block Jacobi preconditioner with incomplete LU factorisation as a subdomain solver. In addition to parallel performance data, we provide new highly-resolved computations of flow around a rotating insect wing and examine its vortex structure and aerodynamic loading.