In this paper, we develop an algorithm to simulate blood flows in aneurysmal arteries and focus on the construction of robust and efficient multilevel preconditioners to speed up the convergence of both linear and nonlinear solvers. The work is motivated by the observation that in the local aneurysmal region, the flow is often quite complicated with one or more vortices, but in the healthy section of the artery, the principal component of blood flows along the centerline of the artery. Based on this observation, we introduce a novel two-level additive Schwarz method with a mixed-dimensional coarse preconditioner. The key components of the preconditioner include (1) a three-dimensional coarse preconditioner covering the aneurysm; (2) a one-dimensional coarse preconditioner covering the central line of the healthy section of the artery; (3) a collection of three-dimensional overlapping subdomain preconditioners covering the fine meshes of the entire artery; (4) extension/restriction operators constructed by radial basis functions. The blood flow is modeled by the unsteady incompressible Navier–Stokes equations with resistance outflow boundary conditions discretized by a stabilized finite element method on fully unstructured meshes and the second-order backward differentiation formula in time. The resulting large nonlinear algebraic systems are solved by a Newton-Krylov algorithm accelerated by the new preconditioner in two ways: (1) the initial guess of Newton is obtained by solving a linear system defined by the coarse preconditioner; (2) the Krylov solver of the Jacobian system is preconditioned by the new preconditioner. Numerical experiments indicate that the proposed preconditioner is highly effective and robust for complex flows in a patient-specific artery with aneurysm, and it significantly reduces the numbers of linear and nonlinear iterations.