Abstract: Camellia oleifera is one of the most important cash crops worldwide. Shelling equipment has been widely used to harvest fresh Camellia oleifera. An accurate discrete element model is often required to characterize the multi-layer mechanical response of the fruit interior during shelling. However, single-layer discrete element models cannot capture the mechanical differences between the shell and the seed, due to treating the fruit as a homogeneous entity. Double-layer models can be expected to distinguish between the outer shell and the internal seed body. Since the seed is treated as a uniform structure, the conventional model cannot capture the sequential failure sequence—from outer shell fracture to seed coat rupture and subsequent kernel response—that occurs under compression. In this study, a multi-layer flexible bonding model was established using the discrete element method and the Hertz-Mindlin with Bonding contact model. Camellia oleifera shared a three-layer composite structure, including the outer shell, the seed coat, and the seed kernel. A bonding mechanism for each layer was employed to characterize material fracture behavior using bond breakage between constituent particles. The basic Hertz-Mindlin contact model was applied between adjacent layers to simulate contact compression and extrusion interaction. Heterogeneous modeling provided a physically realistic representation of intra-layer fracture and inter-layer separation, compared with conventional double-layer models. Meanwhile, geometric parameters—including transverse diameter, longitudinal diameter, and shell thickness—and physical properties —including density, Poisson's ratio, and shear modulus—of fresh Camellia oleifera were determined after systematic physical experiments. Contact characteristic parameters, including coefficients of restitution, static friction, and rolling friction for various material pairings, were calibrated to combine experimental measurements with EDEM simulation. The bonding parameters were calibrated for the seed coat–seed coat, kernel–kernel, and shell–shell contacts. A multi-stage experimental protocol included Plackett Burman screening, steepest ascent climbing, and response surface optimization. In the tea seed model, the relative errors between simulated and experimentally measured peak cracking forces in the transverse and longitudinal directions were 3.17% and 7.13%, respectively. In the whole fruit model, the relative errors were 1.54% and 1.83%, respectively, indicating that the calibrated model accurately predicted the critical failure loads. The measured and simulated compression load displacement curves were compared to evaluate model fidelity beyond single-point peak force in the entire loading history. The coefficient of determination (R²) was employed to quantify the overall agreement between the two curves at the elastic deformation, crack initiation, and crack propagation stages. In the whole fruit, the R² values were 0.92 and 0.90, respectively, for transverse and longitudinal compression; In the tea seeds, the R² values were 0.89 and 0.87, respectively. The high R² values confirmed that the simulation matched the peak failure loads to reproduce the full force deformation response with excellent fidelity. Furthermore, macro cracks in the shells consistently appeared along the longitudinal valve lines in both physical experiments and simulations, indicating the correct failure morphology. The progressive cracking behavior was characterized from the fruit shell and the layer by layer mechanical response within the tea seeds—from seed coat rupture to kernel deformation. The findings can provide a robust theoretical basis and support for the simulation driven optimization of Shelling equipment for Camellia oleifera.