This study presents a novel numerical framework for evaluating damage progression in check dams subjected to impacts from boulder-enriched debris flows. The model integrates three core components: a) a depth-averaged two-phase model to simulate debris flow dynamics accounting for fluid-solid interactions; b) a probabilistic module to represent the stochastic nature of boulder impacts, including variations in size, velocity, and frequency; and c) an elastoplastic contact-based damage model to assess structural degradation under repeated boulder impacts. A finite volume method augmented with grid-optimization techniques is utilized to solve the coupled equations efficiently and accurately resolve localized damage. Validation against experimental data, including flume tests and impact trials from the USGS and bridge pier studies, confirms the model's capability to replicate debris flow behavior and impact forces with high fidelity. Numerical case studies highlight the critical influence of boulder size and velocity (linked to debris flow viscosity) on impact forces: debris flows with lower viscosity (dilute) generate higher peak loads than their more viscous counterparts, leading to distinct patterns of dam damage. Additionally, the sequence of boulder impacts significantly affects both the onset and progression of dam failure. The study emphasizes the necessity of incorporating boulder randomness and debris flow rheology in dam design and hazard assessment. Limitations include simplified assumptions regarding boulder kinematics and dam material properties, suggesting the need for further refinement. This work offers a predictive tool for assessing check dam performance under debris flow hazards, supporting the design of resilient infrastructure in mountainous regions.