We introduce a robust and efficient preconditioner for a hybrid Newton–GMRES method for solving the nonlinear systems arising from incompressible Navier–Stokes equations. When the Reynolds number is relatively high, these systems often involve millions of degrees of freedom (DOFs), and the nonlinear systems are difficult to converge, partially due to the strong asymmetry of the system and the saddle‐point structure. In this work, we propose to alleviate these issues by leveraging a multilevel ILU preconditioner called HILUCSI, which is particularly effective for saddle‐point problems and can enable robust and rapid convergence of the inner iterations in Newton–GMRES. We further use Picard iterations with the Oseen systems to hot‐start Newton–GMRES to achieve global convergence, also preconditioned using HILUCSI. To further improve efficiency and robustness, we use the Oseen operators as physics‐based sparsifiers when building preconditioners for Newton iterations and introduce adaptive refactorization and iterative refinement in HILUCSI. We refer to the resulting preconditioned hybrid Newton–GMRES as HILUNG. We demonstrate the effectiveness of HILUNG by solving the standard 2D driven‐cavity problem with Re 5000 and a 3D flow‐over‐cylinder problem with low viscosity. We compare HILUNG with some state‐of‐the‐art customized preconditioners for INS, including two variants of augmented Lagrangian preconditioners and two physics‐based preconditioners, as well as some general‐purpose approximate‐factorization techniques. Our comparison shows that HILUNG is much more robust for solving high‐Re problems and it is also more efficient in both memory and runtime for moderate‐Re problems.