This thesis presents a novel simulation framework for analyzing the complex dynamical environment surrounding small irregular celestial bodies. A Random Gaussian Sphere generation algorithm is developed to create realistic asteroid shape models, overcoming the limitation of relying solely on pre-existing asteroid data. The framework implements a polyhedral gravity field model to compute gravitational potentials and accelerations, supporting both uniform and non-uniform internal density distributions. The study investigates the combined effects of irregular shape, rotation, and density on orbital and surface dynamics. In particular, evaluating the ridge line around the asteroid and locating the equilibrium points that lie on it, permits to analyze the stability of motion near equilibria through the examination of their eigenstructure, and assess of how physical parameters (density, rotation rate) influence the configuration of the equilibria.