This work integrates stochastic theory with a phase-field model for spallation in ductile metals. By assigning four distinct random distributions to the initial yield strength to characterize the random distribution of material defects and employing an explicit dynamic solver, the entire process of spall damage—from gradual evolution to instability and coalescence—was successfully simulated. The simulation results were validated through plate impact experiments and triangular wave loading experiments. These validations revealed the relationship between the heterogeneity of material yield strength and both the spall strength and the number/area of damaged zones. The results indicate a negative correlation between the standard deviation of the initial yield strength and the spall strength, which holds for both single and multiple spall scenarios in ductile metals. For single spallation, regardless of the initial distribution of yield strength, the resulting spall strength follows a normal distribution. For multiple spallation, the number of initially nucleated damaged zones increases linearly with the standard deviation, while the size of these zones follows a Weibull distribution. Under the same initial random distribution, the number of damaged zones evolves over time, showing a trend of initial slow growth, subsequent acceleration until saturation, and a final decline after saturation. This trend corresponds to the typical process of damage evolution involving nucleation and coalescence during spallation.