Using multi-GPU in lattice Boltzmann method (LBM), fully developed turbulent flow in a square duct at the friction Reynolds numbers (Reτ) of 300, 600, 1200 and 1800 are simulated. Through simulation of three-dimensional lid-driven cavity flow under different Reynolds number (Re), the accuracy of lattice Bhatnager-Gross-Krook (LBGK) multi-GPU program is validated. For turbulent flow in a square duct, all mean velocity, secondary flows, root mean square (rms) of pulsating velocity and Reynolds shear stress predicted by LBGK under the lower Reτ agree well with the literature results, which further verified the effectiveness of the LBGK. In addition, fully developed turbulent flow in a square duct with Reτ up to 1800 predicted by LBGK with 600 million grids provides a reference for turbulent flows under high Reτ . Numerical results show that the LBGK model with low accuracy successfully captures turbulent characteristics for flows at high Re by increasing the grid size, indicating the feasibility and practicality of multi-GPU LBM for modeling industrial flows.