Compute unified device architecture (CUDA) and message-passing-interface (MPI) were used to implement the lattice Boltzmann method (LBM) multi-GPUs parallel algorithm on supercomputing system Mole-8.5E. The accuracy and effectiveness of the multi-GPUs parallel LBM algorithm was verified by the three-dimensional lid-driven cavity flow. Using this parallel algorithm, large-scale simulations of fully developed turbulent square duct flows with Reynolds numbers Reτ of 300, 600, and 1200 were carried out. Numerical results showed that when the grid size was less than the viscous sublayer, Δ+<5, the statistical error of the transfer characteristics near the wall was low, and the simulation accuracy met the need of engineering application. Meanwhile, when Reτ was 300 and 600, the simulation results at different grid sizes Δ+ showed that the LBM had a weak grid-dependent in the statistics of turbulence characteristics in the square duct central area, and when Reτ was 600, the averaged errors of the mean streamwise velocity predicted by LBM at Δ+=1.667, 3.750, and 6.250 compared with direct numerical simulation (DNS) were 1.357%, 2.994%, and 4.766%, respectively. The characteristics of turbulent square duct flow were studied when Reτ was 300, 600, and 1200, and the corresponding grid size Δ+ was 0.833, 1.667, and 3.333, respectively. The secondary flows were successfully captured, and the predicted mean streamwise velocity, pulsating root-mean-square (rms) velocity and other trends and results were consistent with the literature, which further verified the reliability of the single relaxation time (SRT) LBM, and the numerical results also provided a reference for understanding the turbulence characteristics of turbulent square duct flow at high Reτ. The simulation of turbulent square duct flow verified the potential of the SRT LBM multi-GPUs parallel algorithm in ultra-large-scale grid computing, and laid the foundation for further realization of the larger-scale numerical simulation required in practical engineering.