Wind power installations have grown steadily in the last decade, with the growth being driven by reduction in the levelized cost of energy (LCOE) enabling it to compete with conventional power generation sources such as coal and natural gas. These improvements stem from an increase in the average hub height and rotor diameter of wind turbines. Larger rotors also experience greater loads during their operation. In order to enable the structure to withstand the high loads, the larger blades are designed to be flexible and have slender profiles. Additional measures, such as prebending and sweeping of the blades also aid in the design process. However, the current engineering models used to model the blade aerodynamics are limited by their inadequacy in accurately capturing the complex flow effects produced as a result of these out-of-plane geometries. Instead, high-fidelity models involving 3D computational fluid dynamics (CFD) that solve the Navier-Stokes equations can be used to accurately capture the flow around wind turbine blades with out-of-plane geometries. By overcoming the deficits inherent in the low-fidelity aerodynamic models, performing design optimizations using high-fidelity flow models that accurately represent the real-world turbine performance would result in realistic optimum designs. But high computational costs and a myriad of numerical issues preclude performing high-fidelity CFD-based design optimizations. Despite the progress made in recent years in the application of CFD-based shape optimization to 3D rotating wind turbine blades, considerable enhancements are required to make it a mainstay in the industry. Some of the areas of improvement include incorporating realistic structural constraints, deep and stable convergence of flow and adjoint variables, turbulence models for realistic load prediction and their adjoint sensitivities, and a complete set of shape and planform design variables covering the whole rotor. The main objective of this project is to make CFD-based aerodynamic shape optimization of wind turbine rotors practically relevant by enhancing the existing optimization methodologies and solving the specific numerical issues related to using CFD for such a purpose. A proposed way to improve the optimization methodology is to implement a variable-fidelity optimization framework. The benefit would be a reduction in the computational cost associated with performing multiple evaluations of the objective function for a large number of shape and planform variables with a high-fidelity model. This would in-turn allow for the inclusion of a larger design space. When supplemented by relevant structural constraints, the shape optimization would result in realistic designs. The shape optimization would be aided by the implementation of an adjoint solver for the 2D version of DTU Wind Energy's in-house CFD code EllipSys, which would then be used as one of the low fidelity models in the variable-fidelity optimization, with the adjoint-enabled 3D version of EllipSys acting as the high-fidelity model. A numerical issue that has been identified in relation to CFD simulation of rotating flow through wind turbines is of achieving convergence of the solution upto machine precision. In order to limit the computational cost, high-fidelity CFD simulations for wind turbines are performed using steady-state RANS equations with a k-epsilon turbulence model, and for uniform, steady in-flow across the rotor. But the flow through a rotating wind turbine is by its very nature unsteady with flow separation occurring for thicker airfoil profiles of the blades near-root and at the root. In this region, the solution is found to either partially converge or not converge at all. This has a negative effect on the accuracy of the gradients computed using the adjoint method, leading to deceptive solutions for the optimization. This would be improved by applying methods that ensure deep and stable convergence of the flow field.