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 102 to 3 × 104. 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 (Φ = 0.25%, 0.5%, 0.7% and 1%), and the porosity of metal foam (ε = 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 (Nuave ). For AR = 3.75, the value of Nuave for all Rayleigh numbers increased by approximately 5.22 times compared to AR = 1.25. Moreover, Nuave 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 (ε = 0.3), the lowest possible nanofluid concentration (Φ = 0), and the highest possible aspect ratio value (AR = 3.75).