Abstract We present full-volume cosmological simulations, using the moving-mesh code arepo to study the coevolution of dust and galaxies. We extend the dust model in arepo to include thermal sputtering of grains and investigate the evolution of the dust mass function, the cosmic distribution of dust beyond the interstellar medium and the dependence of dust-to-stellar mass ratio on galactic properties. The simulated dust mass function is well described by a Schechter fit and lies closest to observations at z = 0. The radial scaling of projected dust surface density out to distances of 10 Mpc around galaxies with magnitudes 17 < i < 21 is similar to that seen in Sloan Digital Sky Survey data, albeit with a lower normalization. At z = 0, the predicted dust density of Ωdust ≈ 1.3 × 10−6 lies in the range of Ωdust values seen in low-redshift observations. We find that the dust-to-stellar mass ratio anticorrelates with stellar mass for galaxies living along the star formation main sequence. Moreover, we estimate the 850 μm number density functions for simulated galaxies and analyse the relation between dust-to-stellar flux and mass ratios at z = 0. At high redshift, our model fails to produce enough dust-rich galaxies, and this tension is not alleviated by adopting a top-heavy initial mass function. We do not capture a decline in Ωdust from z = 2 to 0, which suggests that dust production mechanisms more strongly dependent on star formation may help to produce the observed number of dusty galaxies near the peak of cosmic star formation.