Recent improvements in accuracy and effciency of numerical simulation techniques in the field of engineering have led to an increasing interest in applying high-fidelity models for wind turbine design. Computational Fluid Dynamics (CFD) based on Reynolds-Averaged Navier-Stokes (RANS) equations in a co-rotating reference frame has shown promising results for performing design analysis. However, using high-fidelity techniques for wind turbine blade design optimization is not yet fully understood. Especially, high-fidelity aerodynamic shape optimization for wind turbine blades has not yet been employed in the industry due to its computational ineffciency when large number of design variables are used with traditional techniques such as finite difference derivatives and gradient-free optimization methods. This dissertation presents an effcient and robust high-fidelity aerodynamic shape optimization methodology for rotating ow problems. The high-fidelity optimization method consists of a multiblock, structured RANS-based CFD simulation tool, a discrete adjoint method, a shape parametrization method, a mesh perturbation technique and a gradient-based optimizer. Steady-state solutions are obtained by the RANS-based CFD analysis method based on a co-rotating reference frame. The turbulence model is a segregated one-equation SA model. Total derivatives of the ow solution and constraint(s) are computed by a discrete adjoint method. For reducing the computational cost of computing partial derivatives that are required in the discrete adjoint method, forward automatic differentiation is used. Once the total derivatives and ow solution are computed, the gradientbased optimizer based on sequential quadratic programming computes a better design by using an augmented Lagrangian formulation with quasi-Newton approximations for the Hessian. The change in design variables obtained by the optimizer is parametrized with a Free Form Deformation (FFD) volume approach. After performing surface perturbations, the mesh is deformed using a hybrid mesh deformation scheme, that combines an algebraic and linear elasticity method. The linear elasticity method based on finite elements is used for large perturbations, while the algebraic method attenuates small perturbations. When the optimality condition is satisfied, the iterative procedure ends with the optimal design. Verification and validation of the developed codes are performed using the NREL VI wind turbine. The RANS-based CFD solver is validated by comparing numerical results with NREL VI sequence S experimental data. The solver resolves attached ow conditions accurately, while separated ow conditions lead to inaccurate ow solutions due to insufficient transition and turbulence modeling. DES can resolve the inaccuracy associated with separated ow conditions. Total derivatives of the adjoint method are verified by comparing derivatives of the complex and finite difference method. The quality of (perturbed) meshes are verified by using mesh quality metrics. Correct shape parametrization is assured after careful examining the direction and magnitude of the deformations. Since the NREL VI wind turbine blade rotates at a constant angular velocity, the power generation is considered to be only dependent on torque. Therefore, the objective of the optimization is maximizing the torque coefficient with shape, twist, and pitch design variables. Thickness constraints between 15% and 50% are added for representing a wing box. The thickness constraints impose thicknesses of the blade to increase only up to 300% of the original thickness. No reduction in thickness with respect to the baseline design at that region of the blade is possible in order to fit the original wing box. For future research, the objective and constraint function(s) can be easily adapted to more realistic and modern rotating ow problems. From the aerodynamic shape optimization of the NREL VI blade, an increase of at least 22.4% in torque is achieved. The airfoil shapes tends to become more cambered and less thick. The nose of the airfoil is more aligned to the in ow. At root section of the wind turbine blade, the trailing edge of the airfoil acts as a ap in order to obtain higher loads at low relative velocity. Three different mesh refinements are employed for optimization. The first mesh is a coarse mesh that is used for verification purposes of the optimization procedure. The second mesh is employed for obtaining accurate aerodynamic shape optimization results. The final design variable values of the medium refined mesh are projected on the most refined mesh, because computational costs would be too high for performing optimization. An increase of 29.1% in torque is achieved, indicating that the increase in optimized torque for more refined meshes will be higher when using coarser meshes. Since wind turbines are operating in a range of wind speeds, multipoint optimization from cut-in to rated wind speed is performed. Similar results as in single-point optimization are achieved. An increase of 22.2% in Annual Energy Production AEP is obtained. The adjoint method and high-fidelity aerodynamic shape optimization methodology allow designers to examine accurately the trade-off between various design variables at the early stage of the design process. For future research purposes, aerostructural and aeroelasticity optimization can be employed with the same framework