A pore-scale numerical method for simulating multiphase granular media is described, in which the solid skeleton is represented by a dense random packing of poly-disperse spherical particles and the pore space is discretized as a network of pores. The underlying idea is to build an efficient method for analyzing the hydro-mechanical behavior of such materials. We design the model specifically for simulating the quasi-static primary drainage, using a criterion for the displacements of the invading fluid and expressions of the forces exerted on the solid grains by fluid’s pressure and surface tension. The model can reproduce some features of drainage, including fluid phase entrapment and capillary fingering. We illustrate the method by simulated primary drainage in oedometer conditions, highlighting the deformation induced by partial saturation.