<p>Space missions to small bodies like asteroids, comets, and moons rely on physics-based simulations to test guidance and control systems. However, accurately modeling their gravitational fields is challenging due to their highly irregular shapes and limited knowledge of their internal structures, complicating orbit planning and landing maneuvers. This study presents a new approach to model realistic density distributions based on Voxel-shaped mass concentrations. We apply body-specific constraints to three-dimensional Perlin noise, supplemented with normalization and segmentation techniques. Additionally, various structural elements can be incorporated into the density distribution. These include centralized and decentralized shells of different thicknesses and densities, as well as anomalies of varying sizes and shapes. Normalization techniques ensure the body’s total mass conservation. We validate our method by calculating the gravitation of a cube and sphere with constant density and comparing it with its analytical solution. We further compare our method with other mascon approaches and the polyhedral method at different Voxel resolutions and conduct additional performance evaluations of our method using test scenarios with focus on geophysical parameters such as the moments of inertia tensor and the gravity field’s spherical harmonics expansion. Our results demonstrate the method’s ability to account for realistic density distributions and to accurately compute the corresponding gravitational fields and geophysical properties.</p>