The concern of this manuscript is stabilized finite element simulations of two-dimensional incompressible and viscous magnetohydrodynamic (MHD) duct flows governed by a coupled system of partial differential equations of the convection–diffusion type. In such flows, the high values of the Hartmann numbers (Ha) indicate that the convection process dominates the flow field. As is typical trouble encountered in convection-dominated flow simulations, classical discretization methods fail to function properly for MHD flows at high Hartmann numbers, yielding several numerical instabilities. In this study, the core computational tool we employ in order to overcome such challenges is the streamline-upwind/Petrov–Galerkin (SUPG) formulation. Beyond that, since the SUPG-stabilized formulation requires additional treatment where solutions experience rapid changes, we also enhance the SUPG-stabilized finite element formulation with a shock-capturing operator for achieving better solution profiles around sharp layers. The formulations are derived for both steady-state and time-dependent models. The proposed method, the so-called SUPG-YZβ formulation, techniques used, and in-house-developed finite element solvers are evaluated on a comprehensive set of numerical experiments. The test computations reveal that the SUPG-YZβ formulation yields quite good solution profiles near strong gradients without any significant numerical instabilities, even for quite challenging values of Ha, whereas the SUPG usually fails alone. Furthermore, by using only linear elements on relatively coarser meshes, the proposed formulation achieves this without the need for any adaptive mesh strategy, in contrast to most previously reported studies.