The reduction in particle size of raw materials using grinding mills is an energy and cost intensive task. Optimization of grinding processes is not trivial as obtaining experimental information is extremely difficult due to the harsh environment. Thus, computational modeling is the most feasible option for obtaining information on the dynamics of the media. However, the computational cost of modeling each particle is high, resulting in the shape of the media being approximated by simple shapes, and in most cases, a reduction in the size of the mill. Even with these simplifications typical simulations take many weeks to months to complete making it infeasible for design prototyping and process optimization. In the last decade, the Graphical Processor Unit (GPU) has enabled large scale simulations of tens of millions of spheres in ball mills using the Blaze-DEM GPU code. Recently, this code was expanded to provide detailed contact detection for polyhedra using the volume-overlap method which is the most accurate approach amongst commercial and academic codes. In this study we first validated the code against experimental results for spherical and cube particle systems in a lab-scale ball mill. Thereafter, we performed a number of ball mill simulations with four additional polyhedral particle systems (truncated tetrahedra, Biluna, elongated hexagonal prisms and a mixed polyhedral particle systems) to study the effect of particle shape. This allows for a first investigation into the roles of particle angularity and aspect ratio on power draw, normal and shear power dissipation between particles, particles and lifters and, particles and the shell. We also show qualitative differences in charge profiles and force chain networks between the various particle systems.