In order to realize the aeroelastic coupling analysis of flexible blade of horizontal axis wind turbine with fewer degrees of freedom, a so-called superelement approach was used to approximate the flexible blade with a number of superelements consisting of a series of rigid bodies connected by springs and dampers. A set of dynamic differential-algebraical equations of blade multibody system was obtained by the R-W method based on the theory of computational dynamics of multibody system. The numerical methods of integration of the set of equations and analysis method for dynamic characteristics of system by impulse response spectrum were presented in this paper. The dynamic response and natural frequencies of a 1.5MW wind turbine blade were analyzed and compared the result with those calculated by ANSYS and ADAMS. The results proved the effectiveness and correctness of the modeling method and the program. The results show that superelement modeling method can describe the elastic deformation of blade and the coupling between dynamic loads and blade deformation well, and the method is a feasible analysis method for the aeroelastic coupling and stability analysis of a wind turbine system.