In this paper, we investigate possible improvements that can be made to the bundle of capillaries model in order to better represent the flow of yield stress fluids through porous media. This was examined by performing extensive and progressive numerical simulations and by introducing the non-circularity of channels' cross section and/or its variability along the channels' axis. It is shown that if only the non-circularity of channels' cross section is taken into account, a moderate influence is observed on both critical pressure gradient for the flow onset and the flow rate/pressure gradient relationship. However, the axial variation in capillaries' cross section has proved to be more impacting the computed flow rate/pressure gradient data. We show hence that when available pore throat and pore body size distributions are used to construct the bundle of axially varying capillaries, the obtained data do fit well experimental results corresponding to the flow of a Bingham-like fluid through a bed of randomly packed mono-sized spheres.