This paper assesses the feasibility of performing topology optimization of laminar incompressible flows governed by the steady-state Navier–Stokes equations using anisotropic mesh adaptation to achieve a high-fidelity description of all fluid–solid interfaces. The present implementation combines an immersed volume method solving stabilized finite element formulations cast in the variational multiscale (VMS) framework and level-set representations of the fluid–solid interfaces, which are used as an a posteriori anisotropic error estimator to minimize interpolation errors under the constraint of a prescribed number of nodes in the mesh. Numerical results obtained for several two-dimensional problems of power dissipation minimization show that the optimal designs are mesh-independent (although the convergence rate does decreases as the number of nodes increases), agree well with reference results from the literature, and provide superior accuracy over prior studies solved on isotropic meshes (fixed or adaptively refined).