A virtual work form of flexible multibody dynamic formulation with rotary inertia has been derived. For the analysis of large flexible multibody systems, deformation modal coordinates have been employed to represent coupled motion between gross and vibrational motion. For the efficient evaluation of the entries in the mass matrix, a flexible body has been treated as a collection of mass points. The rotary inertia was generated from the consistent mass matrix in a finite element model. Deformation mode shapes were obtained from finite element analysis. Bending and twisting vibration analyses of a cantilever have been carried out to see rotary inertia effects. A space flexible robot simulation has been also carried out to show effectiveness of the proposed formulation. This formulation is effective to the model that consists of beam, plate, or shell element that contains rotational degree of freedom at the nodal point. It is also effective to the flexible body model to which a large lumped rotary inertia is attached.