Abstract Two-dimensional numerical simulations were conducted for the natural convection phenomena in a cold square enclosure containing four hot circular cylinders with Rayleigh number in the range of 103 ≤ Ra ≤ 106 and a Prandtl number of Pr = 0.7. The immersed boundary method (IBM) was used to capture the virtual wall boundary of the cylinders based on the finite volume method (FVM). In the present study, the variation in the vertical position of the cylinders was compared with those in the diagonal and horizontal positions. The change in the flow and thermal structures and corresponding heat transfer characteristics were investigated in regard to the transition of the flow regime from the steady state to the unsteady state at a relatively high Rayleigh number. In addition, the effects of the cylinder positions of a rectangular array in the square enclosure on the heat transfer characteristics are highlighted.