Abstract: Cement stabilized materials (CSM) are widely used in pavement base layers, where fatigue damage is inevitable throughout their service life. Due to their significant heterogeneity as multiphase composite materials, a multi-scale approach is essential for studying their fatigue damage. This paper aims to propose a combined approach of simulations and measurements. This approach can characterize the multiphase and heterogeneity properties of CSM and reveal their fatigue damage rules. Firstly, the strength and fatigue performance of CSM were tested, leading to the development of a macroscopic modulus fatigue damage model. Secondly, a pre- and post-fatigue test X-ray computed tomography method with maintained loading was developed to capture the internal meso-structure of CSM. The trainable Weka segmentation was used to provide an accurate meso-structure of CSM for discrete element model (DEM). Thirdly, microscopic testing results were utilized to calibrate the contact parameters of the DEM. The virtually generated aggregate methods for DEM were proposed to enrich the specimens. Finally, virtual fatigue tests were conducted to investigate the fatigue damage rules and to extend the macroscopic modulus fatigue damage model. The results revealed that the fatigue damage rules of CSM accumulate nonlinearly. From a macroscopic perspective, the decay in modulus follows an S-curve across three stages. From a mesoscopic perspective, the average radius coefficient of DEM bonded contacts decreases at an accelerating decay rate. The method proposed in this study reveals the fatigue damage rules under varying stress ratios and cement contents, and develops a simulation based fatigue life prediction equation of CSM. This study offers a reliable numerical technique for modeling and analyzing the fatigue damage rules of composite materials.