Based on the multi-GPU lattice Boltzmann method with the half-way bounce-back scheme, fully developed turbulent duct flows at the friction Reynolds numbers Reτ of 300, 600, 1,200, 1,500, 1,800, and 2,000 were simulated. The parallel performance of multi-GPU lattice Boltzmann simulations is up to 300.162 GLUPS using 1.57 billion grids with 384 GPUs. The simulated friction factor f was consistent with other DNS and experiment results, as well as the Karman–Prandtl theoretical friction law, which verified a sufficient grid resolution Δ+≤3.3, and the LBGK model is stable for Δ+≤5 at high Reynolds numbers. The secondary flows were successfully captured, and turbulence statistics of root-mean-square (r.m.s.) velocity and Reynolds stress were analyzed. The two-point velocity correlation functions and turbulent energy spectra at different positions showed that secondary flows in the near-corner region changed spatial turbulence distribution. Multi-GPU lattice Boltzmann simulations with large grid scales can deal with turbulent square duct flows at high Reynolds numbers and show promise for high-fidelity and scale-resolving fluid dynamics. •Turbulent duct flows are simulated with the multi-GPU lattice Boltzmann method.•Simulated friction factor f is close to the value of the Prandtl–Karman friction law.•The parallel performance of 1.57 billion grids with 384 GPUs is 300.162 GLUPS.•Turbulence and secondary flow in square duct flow are analyzed up to Reτ= 2,000.