A numerical model was used to simulate natural convection heat transfer in a two-dimensional square enclosure filled with nanofluid and saturated with porous media. The enclosure contained three heated tubes with Rayleigh numbers (Ra) ranging from 10(2) to 3 x 10(4). The governing equations with Boussinesq approximation for natural convection were solved iteratively using the finite volume technique through ANSYS FLUENT-CFD commercial package. The Darcy-Forchheimer-Brinkman and Local Thermal Equilibrium models were applied for water and nanofluid flow across the porous zone. The numerical analysis was executed systematically for significant parameters that have a substantial impact on natural convection heat transfer, namely, aspect ratio (AR = 1.25 - 3.75), nanoparticle volume fraction (& phi; = 0.25%, 0.5%, 0.7% and 1%), and the porosity of metal foam (& epsilon; = 0.3, 0.4, 0.5, 0.6 and 0.7). The aspect ratio was found to have a considerable impact on the averaged Nusselt number (Nu(ave)). For AR = 3.75, the value of Nu(ave) for all Rayleigh numbers increased by approximately 5.22 times compared to AR = 1.25. Moreover, Nu(ave) remained approximately constant as the concentration of alumina-water nanofluid increased from 0.25% to 1%. The novelty of this study lies in the multi-objective optimal design approach adopted by adopting DX-12 with the outcomes of ANSYS FLUENT-CFD. The averaged Nusselt number and convection heat transfer rate reached their optimal values by nearly 13 and 9 times, respectively, at the lowest possible porosity (& epsilon; = 0.3), the lowest possible nanofluid concentration (& phi; = 0), and the highest possible aspect ratio value (AR = 3.75).