Sand production in oil wells is often predicted using continuum fluid-coupled models. However, a continuum approach cannot capture important features of the sanding problem, such as erosion and localised failure. This shortcoming of continuum-based analyses can be overcome using the particulate discrete element method (DEM). However, these models, apart from issues of computational cost, have the disadvantage of being difficult to calibrate. One way forward is to calibrate DEM models to capture the response observed in continuum models, where the material parameters can be selected with greater confidence. Adopting this philosophy here, a 3D numerical model based on DEM coupled with Computational Fluid Dynamics was built to simulate sand production around perforations. In the first instance, the basic DEM model (i.e. a dry case) is calibrated against a well-known poro-elastoplastic analytical solution by Risnes et al. (1982). Subsequently, a range of hydrostatic scenarios involving different levels of pore pressure and effective stress are considered. The numerical model shows an asymmetry of the eroded zone that is related to initial microscale inhomogeneity. The stress peak of the analytical solution at the elastic-plastic interface is smoothed because of that asymmetry. The presence of hydrostatic fluid decreases the plastic region and reduces the amount of sand produced. This is not due to changes in effective stress but rather by the particle-scale stabilizing effect of the fluid drag.