In this work, we develop a novel high-order discontinuous Galerkin (DG) method for solving the incompressible Navier-Stokes equations with variable density. The incompressibility constraint at cell interfaces is relaxed by an artificial compressibility term. Then, since the hyperbolic nature of the governing equations is recovered, the simple and robust Harten-Lax-van Leer (HLL) flux is applied to discrete the inviscid term of the variable density incompressible Navier-Stokes equations. The viscous term is discretized by the direct DG (DDG) method, the construction of which was initially inspired by the weak solution of a scalar diffusion equation. In addition, in order to eliminate the spurious oscillations around sharp density gradients, a local slope limiting operator is also applied during the highly stratified flow simulations. The convergence property and performance of the present high-order DDG method are well demonstrated by several benchmark and challenging numerical test cases. Due to its advantages of simplicity and robustness in implementation, the present method offers an effective approach for simulating the variable density incompressible flows.