The modeling and approximation of incompressible flows with variable density are important for a large range of applications in engineering and geophysics. Our main goal here is to develop and analyze novel numerical methods that are suitable for high order finite element and spectral methods. First, we introduce a semi-implicit scheme based on projection methods that uses the momentum, equal to the product of the density with the velocity, as primary unknown. We establish the conditional stability of the scheme and a priori error estimate in time. A fully explicit version of the scheme, that uses time independent stiffness matrices, is introduced and investigated numerically over various setups. Second, we introduce a numerical method based on artificial compressibility techniques that is shown to be more robust than the projection-based method while preserving similar stability and convergence properties. Applications to magnetohydrodynamics instabilities in industrial setups such as aluminum production cells will be presented shortly to show the robustness of the methods.