The nonlinear dynamic problems of three dimensional flexible multibody systems are investigated. The elastic deformation fields of flexible space beams are decomposed into axial deformation and bending deformation, and described by each exact vibration modes in the body coordinate systems. The constrainted nonlinear dynamic equations are derived by using Lagrange multiplier method. A numerical procedure for solving the resulting differential algebraic equations is presented based on Newmark direct integration method combined with the modified Newton-Raphson iterative method. Numerical results verify the effectiveness of the proposed method.