Computational fluid dynamics (CFD) is increasingly used to analyze wind turbines, and the next logical step is to develop CFD‐based optimization to enable further gains in performance and reduce model uncertainties. We present an aerodynamic shape optimization framework consisting of a Reynolds‐averaged Navier Stokes solver coupled to a numerical optimization algorithm, a geometry modeler, and a mesh perturbation algorithm. To efficiently handle the large number of design variables, we use a gradient‐based optimization technique together with an adjoint method for computing the gradients of the torque coefficient with respect to the design variables. To demonstrate the effectiveness of the proposed approach, we maximize the torque of the NREL VI wind turbine blade with respect to pitch, twist, and airfoil shape design variables while constraining the blade thickness. We present a series of optimization cases with increasing number of variables, both for a single wind speed and for multiple wind speeds. For the optimization at a single wind speed performed with respect to all the design variables (1 pitch, 11 twist, and 240 airfoil shape variables), the torque coefficient increased by 22.4% relative to the NREL VI design. For the multiple‐speed optimization, the torque increased by an average of 22.1%. Depending on the CFD mesh size and number of design variables, the optimization time ranges from 2 to 24h when using 256 cores, which means that wind turbine designers can use this process routinely. Copyright © 2016 John Wiley & Sons, Ltd.