To simulate turbulent convection at high Rayleigh number (Ra), we propose a new thermal lattice-BGK (LBGK) model based on large eddy simulation (LES). Two-dimensional numerical simulations of natural convection with internal heat generation in a square cavity were performed at Ra from 10(6) to 10(13) with Prandtl numbers (Pr) at 0.25 and 0.60. Simulation results indicate that our model is fit to simulate high Ra flow for its better numerical stability. At Ra = 10(13), a global turbulent has occurred. With a further increase in Ra, the flow will arrive in a fully turbulence regime. The Nusselt-Rayleigh relationship is also discussed. (c) 2006 Elsevier Ltd. All rights reserved.