Abstract: During the vibration process of wind turbine blades, aerodynamic damping exhibits pronounced temporal and spatial variations. However, due to the limitations of conventional computational methods and the demand for engineering efficiency, its three-dimensional time-varying characteristics are often neglected in practice, and simplified empirical values are instead employed. In this study, a nonlinear dynamic model of the blade is developed by coupling an improved Blade Element Momentum (BEM) method with the Geometrically Exact Beam Theory (GEBT) through Legendre spectral finite elements, and the nonlinear aeroelastic equations are subsequently linearized based on the variational principle. Building on this framework, the Harmonic Balance Method (HBM) is embedded into the linearized model to achieve periodic reduced-order solutions of the system at each harmonic frequency. This enables the decoupled solution of the vibration and flow-field problems, and establishes a novel and efficient approach for computing the three-dimensional time-varying aerodynamic damping of wind turbine blades. The method is validated against aeroelastic wind tunnel experimental data, with the sources of numerical errors analyzed to confirm its accuracy and reliability. Furthermore, the IEA 15 MW reference blade is investigated under typical tip-speed ratio conditions to systematically examine the three-dimensional time-varying characteristics of aerodynamic damping in rotating flexible blades, and the underlying mechanism is elucidated from the perspective of blade surface flow features. The results demonstrate that the proposed method not only efficiently captures the spatiotemporal evolution of aerodynamic damping, but also reveals its complex coupling effects in terms of the dynamic distribution of positive and negative aerodynamic work along the span, the evolution of local vortex structures, and the alternating aerodynamic loads on the pressure and suction sides. Overall, the BEM-GEBT-HBM framework provides significant advantages in both accuracy and computational efficiency, offering an effective tool for the aeroelastic design, optimization, and safety assessment of ultra-long flexible wind turbine blades.