The coupled models of LBM(Lattice Boltzmann Method) and RANS (Reynolds-Averaged Navier-Stokes) are more practical for the transient simulation of mixing processes at large spatial and temporal scales such as crude oil mixing in large-diameter storage tanks. To keep the efficiency of parallel computation of LBM, the RANS model should also be explicitly solved; whereas to keep the numerical stability the implicit method should be better for RANS model. This article explores the numerical stability of explicit methods in 2D cases on one hand, and on the other hand how to accelerate the computation of the coupled model of LBM and an implicitly solved RANS model in 3D cases. To ensure the numerical stability and meanwhile avoid the use of empirical artificial limitations on turbulent quantities in 2D cases, we investigated the impacts of collision models in LBM (LBGK, MRT) and the numerical schemes for convection terms (WENO, TVD) and production terms (FDM, NEQM) in an explicitly solved standard k-epsilon model. The combination of MRT and TVD or MRT and NEQM can be screened out for the 2D simulation of backward-facing step flow even at Re = 10(7). This scheme combination, however, may still not guarantee the numerical stability in 3D cases and hence much finer grids are required, which is not suitable for the simulation of industrial-scale processes. Then we proposed a new method to accelerate the coupled model of LBM with RANS (implicitly solved). When implemented on multiple GPUs, this new method can achieve 13.5-fold acceleration relative to the original coupled model and 40-fold acceleration compared to the traditional CFD simulation based on Finite Volume (FV) method accelerated by multiple CPUs. This study provides the basis for the transient flow simulation of larger spatial and temporal scales in industrial applications with LBM-RANS methods. (C) 2017 The Chemical Industry and Engineering Society of China, and Chemical Industry Press. All rights reserved.