We have developed a new type of fragmentation algorithm that was inspired by a theoretical question raised by A.N. Kolmogorov—and still unanswered after 60 yr—regarding the characteristics of fragment size distributions when the size of the fragments, rβ, exhibits a power-law dependence on the size of the original material, r, with 0 ≤ β ≤ 1. Our fragmentation algorithm uses β and N (which denotes the number of particles produced in the fragmentation) as input parameters and was used for various simulations performed with N values of 2, 3, and 4 and β values from 0 to 1. Simulations with β = 0 resulted in lognormal distributions according to the Kolmogorov–Smirnov goodness-of-fit test at a confidence level of 95%. On the other hand, simulations with fractional values of β > 0 gave highly heterogeneous distributions exhibiting multifractal characteristics. Rényi dimensions (Dq) and Hölder exponents [α(q)] at q = 0 and 1 were defined with coefficients of determination R2 > 0.95 in 78.3% of samples. The resulting α(0), box dimension D0, and entropy dimension D1 = α(1) values spanned the ranges 0.55 to 1.82, 0.52 to 1, and 0.48 to 0.94, respectively, and were thus suggestive of a multifractal nature in the simulated fragment size distributions. The multifractal characteristics of the simulated distributions were consistent with similar analyses performed on actual soil particle size distributions. These results suggest that the new algorithm can be useful for modeling natural fragmentation processes.