To assess the effectiveness of subgrid scale (SGS) models on the prediction results of unsteady loads and turbulent fluctuation of pumpjet propulsors equipped with both front and rear stators, a pumpjet propulsor computational model with attached parts at the model scale is developed using a fully structured mesh, and large eddy simulations are conducted. The computational results of the different SGS models are compared based on five aspects: open water characteristics, turbulence parameters, incoming turbulence spectrum, vortex structure, and fluctuating pressure. Their results are also compared with the experimental values, and the correlation between the internal flow characteristics of the pumpjet propulsor and the turbulent fluctuation is analyzed. According to the results, as regards the prediction of the open water performance of the pumpjet propulsor containing both front and rear stators, the overall trend obtained by the three subgrid models is similar, and the error between the values predicted by the SL model and the experimental ones is the smallest. At the same mesh level, the turbulent fluctuating scale obtained by the SL model is larger than that obtained by the WALE and DSL models, and the turbulent time scale obtained by the DSL model has the smallest fluctuation in the circumferential direction. Among the three SGS models, the turbulent fluctuating scale of the SL model is larger than those of the WALE and DSL models. The SL model exhibits the largest energy dissipation among the three SGS models, followed by the DSL model, while that of the WALE model is the smallest. In the WALE model, the leakage vortex at the top of the blade is the longest, followed by the DSL model, while it is the shortest in the SL model. In the WALE and DSL models, the fluctuating load fluctuates more in the transition region from the middle section to the trailing edge of the blade.