[Objective] In numerical simulations of deep tunnels based on strain gradient theory, the infinitely remote boundary is commonly approximated using a thick-walled cylinder model with a finite outer radius. However, the quantitative influence of the selected outer boundary radius on the computational accuracy of surrounding rock deformation, as well as the existence of a critical radius ratio that ensures computational reliability, remain unclear. This study quantifies the effect of the radius ratio (k) on the computational accuracy of the thick-walled cylinder model for deep tunnels within the framework of Toupin-Mindlin strain gradient theory, with emphasis on characterizing the associated scale effect. [Methods] A plane-strain thick-walled cylinder model based on Toupin-Mindlin strain gradient theory was established. The governing equation was a fourth-order linear homogeneous ordinary differential equation, with boundary conditions involving both Cauchy stresses and higher-order stresses. To solve the high-order boundary value problem, the governing equation was transformed into a system of first-order ordinary differential equations. Using the MATLAB ODE45 solver combined with the quasi-Newton iteration method, an iterative feedback scheme was developed to simultaneously satisfy the stress boundary conditions at the inner and outer boundaries. The validity of the numerical strategy was verified by analyzing the behavior of the model solution as the intrinsic length scale approached zero. By varying k, the strain-gradient numerical solution was compared with the classical elastic analytical solution, and the computational accuracy was quantified using the mean absolute error. [Results] Numerical algorithm showed stable convergence in solving this high-order boundary value problem. As the intrinsic length scale diminished, the numerical solution consistently approached the classical elastic solution, thereby validating the reliability of the proposed methodology. strongly affected the computational accuracy. The numerical error decreased progressively with increasing. When k≥10, the absolute error of the model was effectively controlled to the order of. Parameter sensitivity analysis indicated that the ratio of in-situ stress to support pressure exerted a more pronounced influence on numerical error than the ratio of elastic modulus to the gradient elastic constant. The critical radius ratio criterion of remained valid across various parameter combinations. Based on comprehensive numerical results, an empirical correlation between model error and radius ratio was established. [Conclusion] The findings demonstrate that significantly influences the computational accuracy of strain gradient models for deep tunnels. When k≥10, the computational error induced by finite-domain approximation can be reduced to an acceptable level for engineering applications. The empirical correlation between and provides a quantitative criterion for boundary configuration and practical guidance for determining the minimum computational domain size required to achieve the target error threshold in strain gradient-based numerical simulations of deep tunnels.