Debris flow, a typical non-Newtonian fluid, exhibits rheological behavior significantly influenced by particle size distribution. Traditional rheological models often struggle with applicability and predictive accuracy in complex particle systems. This study proposes a modified Herschel–Bulkley–Papanastasiou (HBP) model, incorporating particle size distribution parameters to dynamically adjust yield stress and shear viscosity, enhancing its accuracy in describing debris flow behavior under varying particle gradations. The model distinguishes the roles of fine and coarse particles: fine particles reduce shear resistance through lubrication effects, while coarse particles enhance yield stress and viscosity via interlocking effects. To validate the model, a series of rheological experiments were conducted on 14 particle gradation conditions. Results showed the modified HBP model achieved fitting coefficients between 0.933 and 0.990, significantly outperforming traditional models and demonstrating superior adaptability across different particle distributions. The model was further integrated into the OpenFOAM framework for three-dimensional simulations of a flume experiment. These simulations considered wall friction and dynamic free surface changes. Comparative analysis with physical experiments revealed the modified HBP model accurately captured debris flow behavior, free surface dynamics, and pressure field distributions under varying channel bed conditions. In summary, the modified HBP model overcomes limitations of traditional models by incorporating particle size distribution parameters, offering a more precise and versatile framework for debris flow rheology. This work provides a robust theoretical and numerical tool for advancing the prediction and mitigation of debris flow in engineering applications.