A new numerical algorithm was developed to provide a fully detailed flow field in liquid metal MHD flow with a relatively large Hartmann number and interaction parameter. The algorithm includes the effects of advection and diffusion, and is capable of predicting momentum and heat transfer in MHD flows. Using this algorithm, an incompressible, viscous, three-dimensional MHD flow in a square duct is investigated at a low magnetic Reynolds number by means of the finite volume method. The velocity and temperature profiles are obtained in the developing region for constant wall temperature. The result shows that large velocities are obtained near the insulating walls parallel to the magnetic field. Also, near the perfectly conducting walls perpendicular to the field, a velocity profile like a Hartmann layer is obtained. In association with the velocity profiles, Nusselt number at the insulating walls (with side layer) is seen to be larger than that at the perfectly conducting walls (with Hartmann layer).