Accurate modeling and efficient solution are crucial for performance analysis and optimization of thermal systems. However, this task is increasingly difficult with the system scale and complexity increase and the coupling of nonlinear system constraints. To deal with this obstacle, most studies introduce simplifications to improve calculation efficiency and stability, while the accuracy of model is sacrificed. This work proposes an efficient and robust algorithm for thermal system simulation based on Generalized Benders Decomposition (GBD) technique. A supercritical carbon dioxide (sCO₂) recompression Brayton cycle is used to present the algorithm in detail. The algorithm is verified by the high agreement with experimental results reported. It is further compared with conventional simulation methods to highlight its advantages on the sCO₂ cycle. The simultaneous equations method fails to deal with the system simulation without simplifications due to the drastic property variation of sCO₂. The sequential modular method (SMM) requires eight layers of iteration, of which four layers are nested, leading to a long computation time. In contrast, the proposed algorithm only has one iteration layer without nested iterations, and hence it consumes less than 10% calculation time compared with SMM. Moreover, the proposed algorithm also owns higher robustness. Finally, the sCO₂ Brayton cycle is optimized to improve its efficiency by integrating the proposed GBD-based simulation method and genetic algorithm. Results show that the system efficiency reaches 48.83% after operation parameter optimization.