We present novel numerical schemes for the incompressible Navier-Stokes equations with variable density and address two critical concerns, i.e., the preservation of lower density bounds and the preservation of energy inequality under the gravitational force. Spatial discretization is performed using upwind discontinuous Galerkin methods for density, and continuous finite element methods are used for velocity and pressure. Additionally, velocity fields are projected onto divergence-free Raviart-Thomas finite element space in the transport equation. This relaxes using divergence-free finite elements and therefore simplifies the implementation. The proposed schemes are tested on benchmark examples to demonstrate the excellent performance in accurately capturing the complex dynamics of these problems. Overall, the proposed numerical schemes exhibit enhanced stability, efficiency, and accuracy, making them suitable for modeling and simulating incompressible flows with variable density in practical applications.