The issue of mesh quality has become a topic of increasing interest among scientists and professionals working in the field of computational fluid dynamics, as it has a major impact on model performance. Here we introduce a mesh generation method that is based on a Voronoi decomposition of the flow domain, resulting in polygonal regions in 2D and polyhedral cells in 3D. The discretised Reynolds equations for arbitrary cells are derived using a fmite volume approach and the required boundary conditions are discussed. Both the grid generation process and a solver for the turbulent flow equations have been implemented in the simulation model RSim-3D. The model has been applied to several test cases and the results have been satisfactory. Polyhedral cells have the potential to delivering more accurate results under certain flow conditions than hexahedral cell types.