Two-dimensional numerical simulations were conducted to investigate the natural convection heat transfer induced by a temperature difference between the hot surfaces of four inner cylinders and the walls of an enclosure. The simulations were done using different Rayleigh numbers (10(3) <= Ra <= 10(6)) and distances between neighboring cylinders (0.3 <= epsilon <= 0.7). The immersed boundary method (IBM) was used to capture the virtual boundary of the four cylinders, which we arranged in a diamond array. Detailed analysis results including the distribution of streamlines, isotherms, and Nusselt numbers are presented. The thermal and flow fields reached steady state in the Rayleigh number range of 10(3) <= Ra <= 10(5), regardless of the variation in epsilon. However, at Ra = 10(6), the numerical solutions show time-dependent characteristics when the distance between neighboring cylinders is 0.4 <= epsilon <= 0.6. (C) 2017 Elsevier Ltd. All rights reserved.