<p>In this paper, we proposed a novel grain-based model based on particle flow code to realistically reproduce the heterogeneous structure of crystalline granite. Then, it is applied to the cyclic loading and unloading simulation. Based on the quantitative analysis of the three-dimensional multilevel force chain network, the evolution of force chain characteristics of crystalline granites with different minimum radii of the grains <i>R</i>_G during cyclic loading and unloading is investigated. Our results demonstrate that specimens with varying <i>R</i>_G exhibit stress–strain curves that form a “hysteresis loop” due to nonideal elasticity deformation. As<i> R</i>_G increases, the proportion of intragranular contacts with higher micro-strength and micro-modulus rises, enabling it to bear more loads and exhibit greater deformation resistance. The microscale slip between particles is also reduced when an intragranular contact fractures. Consequently, both the upper stress threshold and the elastic modulus of the sample increase as<i> R</i>_G increases, while the variation range of strain values decreases. During loading, most cracks primarily propagate in an orientation range orthogonal to the loading direction. As<i> R</i>_G increases, the average value and sum value of whole general force chains increase. The main orientation of high-strength force chains (HF) aligns with the loading direction. With an increase in <i>R</i>_G, the numbers of HF in whole structures and intragranular structures rise, while the number of HF in intergranular structures decreases. As <i>R</i>_G increases, the number of basic elements and contacts within the mineral structure that can jointly bear the load increases, and the formed force chain network can bear a higher level of load. Due to the difference of micro-strength, the bearing capacity of the intragranular structure is greater than that of the intergranular structure.</p>