<p>The divergence constraint of the incompressible fluids usually causes the weak robustness of standard mixed finite element methods. Grad-div stabilization is a popular technique for improving the robustness. In this paper, we theoretically show that for magnetohydrodynamic flows at large Hartmann numbers, grad-div stabilization can improve the well-posedness and robust stability of the continuous problem, and remove the effect of Hartmann number on the finite element discrete errors. Besides, applying the backward Euler method and lagging the nonlinear term, we construct a linear grad-div stabilized finite element algorithm for magnetohydrodynamics flows at low magnetic Reynolds numbers. A complete theoretical analysis of its stability and convergency is provided. Some computational experiments illustrate the validness of our algorithm and its theoretical results and also the benefits of grad-div stabilization.</p>