<p>Direct numerical simulations have been conducted to investigate the evolution process of liquid metal laminar to turbulent flow in a rectangular duct under the influence of a non-uniform magnetic field. The Reynolds number is <i>Re</i> = 6299, and the inlet Hartmann number is <i>Ha</i> = 2900, with the magnetic field strength decreasing along the flow direction. The results indicate that the dynamic reversal of the three-dimensional (3D) Lorentz force direction near the inflection point of the magnetic field dominates the flow reconstruction, driving the wall jet acceleration and forming an M-type velocity distribution. Moreover, the high-speed shear layer of the jet triggers Kelvin-Helmholtz instability, resulting in the generation of secondary vortex structures near the parallel layer in the non-uniform magnetic field region. In the cross-section perpendicular to the flow direction, the secondary flow gradually evolves into a four-vortex structure, while the velocity fluctuations and turbulent kinetic energy reach the peak. Based on the characteristics of the vortex rotation direction near the shear layer, the intrinsic mechanism behind the unique bimodal distribution of the root-mean-square of velocity fluctuations in the parallel layers is revealed. Furthermore, by comparing the evolution of turbulence under different magnetic field gradients, it is revealed that the distributions of shear stress, Reynolds stress, and turbulent kinetic energy exhibit significant parameter dependence. The strong 3D magnetohydrodynamic effects at the magnetic field gradient <i>γ</i> = 0.6 have an immediate impact on the pressure distribution. The transverse Lorentz force <i>LF</i>_<i>z</i> further promotes the fluid to accumulate at the wall, leading to a significant increase in the pressure drop and transverse pressure difference in the flow.</p>