The mechanical response of geotechnical granular materials exhibits multi-scale characteristics. The multi-scale simulations coupling the Finite Element Method (FEM) and Discrete Element Method (DEM) can effectively capture the multi-scale responses while maintaining high computational efficiency. A GPU-parallel FEM-DEM coupling code was developed based on the high-performance DEM software MatDEM, with computational parameters and results analyzed in conjunction with granular pore fractal characteristics. Firstly, multi-fractal theory was employed to investigate the spatial distribution characteristics of pores, identifying key fractal indices. Subsequently, the reliability of the FEM-DEM coupling code was verified through single element tests and biaxial compression tests. Finally, the meso-scopic responses of Representative Volume Elements (RVEs) at different locations were investigated based on biaxial compression tests. Results demonstrate that the pore spatial distribution within RVE exhibits multi-fractal characteristics. When the particle quantity exceeds 400, the self-similarity of pore spatial distribution ensures the stability of stress-strain responses output after homogenization of granular assemblies within RVE. The capacity dimension D₀ and singularity index α₀, which characterize the average information of pore distribution, show linear correlations with RVE volumetric strain. These indices can serve as internal variables reflecting the complexity of granular material spatial characteristics. This study provides an exploration for analyzing the macro-meso mechanical relationships in engineering-scale granular deposits.