Natural convection heat transfer from vertical 99 rod bundle in liquid sodium was numerically analyzed to optimize the thermal-hydraulic design for the bundle geometry with equilateral square and triangle arrays (ESA and ETA). Unsteady laminar three-dimensional fundamental equations of natural convection heat transfer caused by step heat flux were solved numerically until the solution reached a steady state. The PHOENICS code was used to calculate the temperature dependence of the relevant thermo-physical properties. For this work, a 99 test rod with a diameter (D = 7.6 mm), heated length (L = 200 mm), L / d (= 26.32) was used. The surface heat flux (q) of each cylinder was equally applied. It ranged from 210 5 to 710 6 W/m 2 . The modified Rayleigh number, (Ra f,L ) ij and (Ra f,L ) NxNy,S/D , ranged from 6.43010 5 to 4.30910 7 at the liquid temperature of 673.15K. The value of S/D, the ratio of the diameter of the flow channel to the diameter of the rod in the bundle geometry of the vertical 99 rod bundle, ranged from 1.82 to 5 for the bundle geometry with ESA. The spatial distributions of the average Nusselt numbers, (Nu av ) ij and (Nu av,B ) NxNy,S/D , were revealed for the vertical single cylinder of the rod bundle and the vertical rod bundle. To determine the effect of array size, bundle geometry, S/D, (Ra f,L ) ij and (Ra f,L ) NxNy,S/D on heat transfer, the average Nusselt numbers, (Nu av ) ij and (Nu av,B ) NxNy,S/D , of bundle geometries with different values of S/D were calculated. The general correlations for natural convection heat transfer from a vertical N x N y rod bundle with the ESA and ETA including the effects of array size, (Ra f,L ) NxNy,S/D and S/D were derived. The correlations of vertical N x N y rod bundles can describe the theoretical (Nu av,B ) NxNy,S/D values within -4.79% to 6.60% at S/D ranged from 1.82 to 6, and modified Rayleigh numbers ranged from 6.43010 5 to 4.30910 7 .