A consistent and conservative scheme designed by Ni et al. for the simulation of MHD flows with low magnetic Reynolds number has been implemented into a 3D parallel code of HIMAG based on solving the electrical potential equation. The scheme and code are developed on an unstructured collocated mesh, on which velocity (u), pressure (p), and electrical potential ([variant phi]) are located in the cell center, while current fluxes are located on the cell faces. The calculation of current fluxes is performed using a conservative scheme, which is consistent with the discretization scheme for the solution of electrical potential Poisson equation. The Lorentz force is calculated at cell centers based on a conservative formula or a conservation interpolation of the current density. We validate the numerical methods, and the parallel code by simulating 2D fully developed MHD flows with analytical solutions existed and 3D MHD flows with experimental data available. The validation cases are conducted with Hartmann number from 100 to 104 on rectangular grids and/or unstructured hexahedral and prism grids.