This work investigates the implementation of nonlinear model reduction to°exible multibody dynamics. Linear elastic theory will lead to instability issues with rotating beam-like structures, due to the neglecting of the membrane-bending coupling on the beam cross-section. During the past decade, considerable e®orts have been focused on the derivation of geometric nonlinear formulation based on nodal coordinates. In order to reduce the computation cost in°exible multibody dynamics, which is extremely important for complex large system simulations, modal reduction is usually implemented to a rotating°exible structure with geometric nonlinearities retained in the model.In this work, a standard model reduction process based on matrix operation is developed, and the essential geometric sti®ening nonlinearities are retained in the reduced model. The time responses of a tip point on a rotating EulerÀBernoulli blade are calculated based on two nonlinear reduced models. The¯rst reduced model is derived by the standard matrix operation from a full¯nite element model and the second reduced model is obtained via the Galerkin method. The matrix operation model reduction process is validated through the comparison of the simulation results obtained from these two di®erent reduced models.An interesting phenomenon is observed in this work: In the nonlinear model, if only quadratic geometric sti±ng term is retained, the two reduced models converge to the full¯nite element model with only one bending mode and two axial modes. While if both quadratic and cubic geometric sti±ng terms are retained in the nonlinear equation, the modal-based reduced model will not converge to the¯nite element model unless all eigenmodes are retained, that is the reduced model has no degree of freedom reduction at all.