As fusion liquid metal blankets work under strong magnetic fields, the liquid metal would form magnetohydrodynamics (MHD) flow of high Hartmann number in ducts. With the anisotropic suppression effect of magnetic field towards turbulence fluctuation, MHD turbulence becomes quasi-2 dimensional, whose characteristics of heat and mass transfer are different from common turbulence. Present MHD turbulence models can not accurately predict this suppression effect. In this paper, a K-Epsilon RANS turbulence model was developed for incompressible MHD flow at high Hartmann number in fusion liquid metal blankets. According to the characteristic magnetic braking time of isotropic turbulence in MHD duct flows, the effectiveness of the electromagnetic damping process was simulated. Turbulence in MHD rectangular duct flow of high Hartmann number and Reynolds number with strong sidewall jets was chosen to benchmark the model. The MHD turbulent velocity and turbulent kinetic energy were compared between the model and the experiment. The benchmark of MHD duct flow showed that the model could accurately simulate the instability location formed by MHD jets, the distribution of turbulent kinetic energy, and variation of MHD quasi-2-dimensional turbulence strength with magnetic fields.