The depletion of fossil fuel resources and the growing need for sustainable energy solutions have increased interest in vertical axis wind turbines (VAWTs), which offer advantages in urban and variable-wind environments but often exhibit limited performance at low tip speed ratios (TSRs). This study optimizes VAWT aerodynamic behavior across a wide TSR range by varying three geometric parameters: maximum thickness position ($a/b$), relative thickness (m), and pitch angle ($\beta$). A two-dimensional computational fluid dynamics (CFD) framework, combined with the Metamodel of Optimal Prognosis (MOP), was used to build surrogate models, perform sensitivity analyses, and identify optimal profiles through gradient-based optimization of the integrated $C_p$–$\lambda$ curve. The Joukowsky transformation was employed for efficient geometric parameterization while maintaining aerodynamic adaptability. The optimized airfoils consistently outperformed the baseline NACA 0021, yielding up to a 14.4% improvement at $\lambda =2.64$ and an average increase of 10.7% across all evaluated TSRs. Flow-field analysis confirmed reduced separation, smoother pressure gradients, and enhanced torque generation. Overall, the proposed methodology provides a robust and computationally efficient framework for multi-TSR optimization, integrating Joukowsky-based parameterization with surrogate modeling to improve VAWT performance under diverse operating conditions.