<p>This study comprehensively investigates the grinding process of spiral bevel gears and establishes a stochastic multi-grain finite element model using ABAQUS. A single-factor simulation approach is employed to systematically examine the depth-dependent distribution of residual stresses within the subsurface following grinding. Based on these results, a predictive model for residual stress is developed using a three-factor, four-level orthogonal experimental design. Grinding parameters are subsequently optimized through Grey Relational Analysis (GRA) and Particle Swarm Optimization (PSO), with a comparative evaluation of both methods conducted via the entropy weight approach. Key findings indicate that the machined tooth flanks exhibit a residual compressive stress state, with the maximum compressive stress located in the near-subsurface region. This stress gradually transitions to residual tensile stress with increasing depth, eventually approaching zero. Among the grinding parameters, depth of cut exerts the most pronounced influence on residual stress magnitude, followed by grinding speed, whereas feed rate has a comparatively minor effect. The GRA-based optimization scheme achieves superior solution quality compared with the PSO approach. Experimental validation demonstrates strong agreement between measured and simulated residual stress profiles, with overall discrepancies below 16%, thereby confirming the reliability of both the multi-grain simulation model and the predictive framework. These findings provide a robust theoretical foundation for optimizing the grinding process in spiral bevel gear manufacturing.</p>