Seepage failure is a common problem in engineering, and the calculation and analysis of critical hydraulic gradient are of great significance for the safety and protection of engineering. Based on the principle of discrete element method and computational fluid dynamics, the fluid–solid coupled models were established to study the critical hydraulic gradient and particle loss rate of granular soils at seepage failure. The evolution of seepage failure was divided into four stages: seepage development stage, local damage stage, volume expansion stage and overall damage stage. The validity of numerical simulation was demonstrated by comparing the critical hydraulic gradient obtained by numerical simulation and by Terzaghi’s formula. According to the fabric damage and flow velocity variation of the models at seepage failure, the influences of model size and particle size on the critical hydraulic gradient and particle loss rate were analyzed. The results indicate that critical hydraulic gradient and particle loss rate were not sensitive to changes in model size. A wide particle size distribution range resulted in large critical hydraulic gradient and small particle loss rate at seepage failure. The discrete element numerical simulation can not only be used to determine the critical hydraulic gradient of geotechnical and hydraulic engineering, but also offer a visual portrayal of the evolution of seepage failure, serving as an important complement to comprehend the intricate microscopic mechanisms underlying soil seepage failure.