High fidelity large-scale direct numerical simulation of particulate flows is of great value in a variety of industrial applications. It is computationally intensive as it combines time integration, solving nonlinear algebraic equations, and the associated linear systems. The finite element discretization of the coupled system of PDEs on an unstructured grid using an arbitrary Lagrangian–Eulerian moving mesh technique leads to very large nonlinear systems that are linearized by a version of Newton's method. The linear algebraic systems (Jacobians) are sparse, nonsymmetric and indefinite, for which standard linear system solvers based on Krylov subspace methods generally fail to converge without appropriate preconditioners. The failure of Krylov methods in production codes is currently being addressed by reducing the size of the time step. This, however, leads to a very long simulation time, and therefore is not always a viable approach. In this study, we design a hybrid inner–outer iterative scheme for solving these indefinite systems which proves to be both efficient, robust and ideally suited for parallel computing platforms even with appropriate large time steps. Comparisons with Krylov subspace methods show the superiority of our proposed class of nested iterative schemes which is also scalable with respect to mesh size, and insensitive to changes in properties of the fluid–particles system.