This report presents a comprehensive study of dark matter physics within the framework of thermal freeze-out and Boltzmann equation formalism. Beginning with cosmological and observational motivations, we derive the relativistic Liouville operator and systematically construct the Boltzmann equation for thermal relics in an expanding universe. Both analytical and numerical methods are employed to solve the evolution equations for particle yield and relic density. We construct simplified dark matter models, focusing on scalar singlet extensions of the Standard Model, including both real and complex scalar fields interacting via the Higgs portal. For these models, we compute the thermal annihilation cross-sections into Standard Model final states and numerically scan the parameter space using Python to determine regions compatible with observed dark matter abundance. The results highlight the role of resonance effects, particularly around the Higgs pole, and demonstrate the reliability of a Python-based approach in probing viable dark matter parameter spaces . This framework lays a foundation for further extensions involving effective operators, direct detection constraints, or non-thermal dark matter production.