We present a high-fidelity aerodynamic shape optimization framework for the NREL S809 airfoil using Free-Form Deformation (FFD) parameterization coupled with a RANS CFD solver and a discrete adjoint for gradient computation. The design objective is to maximize aerodynamic efficiency (lift-to-drag ratio, CL/CD) under single- and multi-point operating conditions relevant to horizontal-axis wind turbines. Geometry changes are controlled by 40 FFD surface control points, and the Sequential Least Squares Programming (SLSQP) algorithm enforces lift, thickness and volume constraints. We validate our solver against NREL Phase VI wind-tunnel data and perform mesh-convergence and gradient-verification studies. Single-point optimization yields a CL/CD increase of ≈22.5% relative to the baseline at Re = 3.48×10⁶, while multi-point optimizations (2–4 points) produce robust improvements across the operating envelope. We discuss sensitivity to turbulence/transition modeling, show gradient check results (finite-difference vs adjoint), and provide mesh-independence evidence. The optimized shapes reduce maximum thickness and move camber forward, improving lift while reducing drag; however, structural implications of thickness reduction are discussed. The framework and validation illustrate a reproducible path for high-fidelity airfoil optimization targeted at wind-energy applications.