In this study, stability analysis of lattice Boltzmann methods has been investigated for simulation of high Ra no. turbulent thermal flow problems and a stable model has been proposed. The following issues related to simulate these problems using lattice Boltzmann method were investigated: total number of links of lattice model, stability improvement techniques, turbulence modeling, compressibility effects, thermal lattice Boltzmann models and implementation of boundary conditions. In two-dimensional case, a combination model has been suggested and applied on a non-uniform mesh. This model is composed of 1. D2Q9 lattice model, 2. Multi-Relaxation-Time Lattice Boltzmann (MRT-LB) model, 3. Incompressible LB model (Guo based), 4. LES in conjunction with Smagorinsky and mixed scale viscosity sub-grid models, 5. Hybrid thermal lattice Boltzmann model, 6. Mid-grid boundary conditions, and 7. Application of the TLLBM technique on non-uniform mesh. To validate the proposed model, a turbulent RBC flow has been simulated at different Ra numbers up to 10 15 for Prantdl number of 0.71 . Besides, in three-dimensional case, a combination model has been proposed. This model consists of 1. D3Q19 lattice model, 2. Fractional Volumetric Multi-Relaxation-Time Lattice Boltzmann (FV-MRT-LB) model, 3. Incompressible LB model (Guo based), 4. LES in conjunction with Smagorinsky and viscosity mixed scale sub-grid models, 5. Hybrid thermal lattice Boltzmann model, and 6. Mid-grid boundary conditions. To verify the proposed model, a turbulent NC flow is simulated at different Ra numbers up to 10 15 for Prantdl number of 0.71 . The results show that the proposed models produce accurate results at low Ra number flows and are stable at high Ra number flow problems. Keywords: Lattice Boltzmann method; thermal LBM; stability of LBM; LES; Rayleigh number.