Laminar and turbulent natural convection inside concentric spherical shells with isothermal cold and hot boundaries is numerically investigated up to Rayleigh number values Ra <= 10(12) and Pr = 0.71. The study utilizes direct numerical simulation (DNS), large eddy simulation (LES) and Reynolds averaged Navier-Stokes (RANS) approaches for investigation of the laminar, transitional and fully developed turbulent flow regimes, respectively. Three-dimensional flow patterns for slightly supercritical oscillatory flow regime inside the shell, with internal/external diameter ratio equal to D-i/D-o = 0.714 are presented and may be potentially useful for verification of the future linear stability analysis results. Particular attention has been given to the complex, fully three-dimensional unsteady flows occurring in narrow shell geometries characterized by 0.85 <= D-i/D-o <= 0.95. For this geometry considerable deviations in predicted heat flux rate through the shell boundaries are observed when compared with existing heat transfer correlations for the entire range of Ra numbers. The deviations tend to increase for transitional and fully turbulent flows. A new correlation for the heat transfer rate is suggested for laminar and transitional flow regimes. (C) 2013 Elsevier Ltd. All rights reserved.