In this paper, we present the development of a rigid-flexible multibody model which, coupled with an existing aerodynamic model, is used to numerically simulate the non-linear aeroelastic behavior of large horizontal-axis wind turbines. The model is rather general, different configurations could be easily simulated though it is primarily intended to be used as a research tool to investigate influences of different dynamic aspects. It includes: i) a supporting tower; ii) a nacelle which contains the electricity generator, the power electronics and the control systems; iii) a hub, where the blades are fixed, connected to the generator rotating shaft; and, iv) three blades which extract energy from the airstream. The blades are considered flexible, and their equations of motion are discretized in space domain by using beam finite elements capable of taking into account the non-linearities coming from the kinematics. The tower is also considered flexible, but its equations of motion are discretized by using the method of assumed-modes. The nacelle and hub are considered rigid, and their equations of motion take into account the effects of the kinematic non-linearities. Due to the system complexity, the tower, nacelle and hub are modeled as a single kinematic chain and each blade is modeled separately. Constraint equations are used to connect the blades to the hub. The resulting governing equations are differential-algebraic, and these are numerically and interactively solved in the time domain by using a fourth order predictor-corrector scheme. The results help to understand the wind speed influence on: i) the rotor angular speed; ii) the after-forward and side-to-side displacements of the tower; and, iii) the flap- and edge-wise displacements of the blades. Copyright (C) 2011, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.