Accurately assessing the erodibility of geomaterials is of great significance for the design of earthen structures and the prevention of the associated failure induced by seepage force. Recently, the un-resolved Computational Fluid Dynamics–Discrete Element Method (CFD–DEM) has been widely used to investigate internal erosion. However, due to the use of wall boundary and the fact that the fixed CFD domain cannot be changed with the soil sample’s volume contraction during the erosion test, a larger porosity at the boundary of the CFD domain is commonly formed, resulting in sidewall preferential flow (i.e., relatively more fine particles migrate along the boundary of the DEM domain) and thereby overestimating the soil erodibility. In this study, a new method based on particle boundary is developed to tackle this problem. The newly proposed particle boundary can prevent its particles from erosion via inter-particle bonding and transfer stress from servo walls to the simulated sample. An optimal particle boundary thickness is determined by considering sample contraction and computational efficiency. The performance of the new method was compared with the conventional method and also verified using experimental results. The results show that the newly proposed method has significantly improved the uniformity of fluid velocity distribution. Furthermore, the cumulative eroded mass of fine particles in the new model is approximately 15% lower than in the conventional model. It is convincingly demonstrated that the new method can simulate internal erosion better and give a more accurate assessment of geomaterial erodibility.