The energy efficiency of recovery boiler power plants is largely influenced by the heat transfer to the superheaters. In the design process of such very large-scale applications, one of the key challenges is the a priori geometry optimization by robust numerical approaches. The main objectives of this work are to demonstrate a numerical optimization framework and to optimize the geometry of the superheater region to enhance the heat transfer. The framework is implemented as a surrogate-based optimization method, which combines simulated annealing, local polynomial regression, and computational fluid dynamics. The novelty of this work consists of the following: 1) The optimization framework is designed and introduced. 2) The connection between the geometry and heat transfer is quantified by formulating the optimal design curve. 3) The optimal design for a typical, existing recovery boiler is identified. The results indicate that the uniformity of the flow field is improved, and the heat transfer rate is increased by 5%. 4) The results indicate the importance of minimizing the separation vortices through the geometrical design with a strong linkage to the overall heat transfer rate. 5) The potential of optimization methods in this very large-scale energy production application is demonstrated for the first time. (C) 2018 Elsevier Ltd. All rights reserved.