This paper presents a meso-scale numerical modelling method to investigate the stability of tunnels in sandy cobble stratum. Using Monte Carlo simulation, Python routine is implemented to randomly generate and filter the location information of blocks including center coordinates, sizes and orientation angles, then extract and import them into CAD modelling software in batches to mesh. The problem of overlapping detection for elliptical and ellipsoidal blocks is solved based on affine transformation, achieving higher efficiency and applicability than the existing methods. Then, using the proposed numerical modelling method, convergence-confinement analysis of circular tunnels with different cover depths is performed in plane strain condition. Meanwhile, face stability analysis of shield tunnels is carried out in three-dimensional condition. Results show that sandy cobble stratum may not be in accordance with the scale-independent soil-rock mixture. Blocks with small sizes are geotechnically significant and contribute much to the overall strength. When the cover depth ratio reaches a certain level, limit support pressure remains basically unchanged due to the strong arching effect above tunnel crown, while permissible support pressure increases approximately linearly to control deformation. Compared with the multi-scale model, hardening soil model considering the overall strength is suitable for the stability analysis of tunnels in sandy cobble stratum at the macro-scale with high precision, while the macro homogeneous model ignoring the presence of blocks provides a conservative (i.e. uneconomical) prediction for the design of tunnel stability.