This thesis presents a novel compaction simulation method, where a multibody dynamics model of a vibratory roller is coupled to the discrete element method (DEM) for unbound granular materials. A multibody dynamics solver with constraints is developed, and an analysis of the roller mechanical system is performed to construct a rigid body model of the machine. Further, a method for bed initialisation is developed, and a uniaxial strain test is used to calibrate the multisphere DEM material parameters. Then, compaction simulation is made possible by implementation of a coupling algorithm, that runs the DEM and multibody solvers simultaneously at different timesteps. Such simulations are compared to full-scale experiments and compaction theory. The machine response to beds of varying stiffness agrees with experience from compaction practice, and characteristic behaviour, such as double jumps, is observed. The stresses in the bed agree with experiments if the particle Young’s modulus is kept low. However, the roller penetration of the bed is higher than in experiments due to insufficient shear resistance in the DEM model. At the same time, no increase in bulk density is achieved. Further analysis shows that the lack of shear resistance is likely related to the multisphere model of the particle geometry, and the lack of compaction may be due to a particle size distribution that is too narrow. On the other hand, such simplifications are necessary, because large computational costs impose limits on the particle size distribution, particle discretisation, and domain size. The function of the solver coupling and machine model is verified, but in order to enhance agreement with full-scale experiments, improvements are needed in terms of both the particle modelling and computational performance capacity.