[en] A finite element analysis of the fully developed laminar natural convection in vertical rod bundles is presented. A Galerkin formulation is used with piecewise linear interpolation polynomials as shape functions and triangular elements. No restrictions are required on radial temperature boundary conditions. The global stiffness matrix is skyline stored. Natural convection in closed hexagonal 7-rod and 91-rod bundles with triangular arrangement is computed for heat fluxes prescribed at the rod surfaces and temperatures specified at the bundle wall. The results show that a coarse mesh is adequate to correctly represent the flow for a large range of Rayleigh numbers. New physical aspects of the natural convection process in finite closed rod bundles are available. This concerns mainly three points. First, ascendant velocities exist in the peripheral region for high Rayleigh numbers. Second, the maximum temperature moves from the center to the periphery of the bundle as the Rayleigh number increases. Third, the use of the half bundle width as a characteristic length is preferable to the hydraulic diameter