SUMMARYDi erent ÿnite-element-based strategies used to represent the components' exibility in multibody systems lead to various sets of co-ordinates. For systems in which the bodies only experience small elastic deformations it is common to use mode component synthesis to reduce the number of generalized elastic co-ordinates and, consequently, the equations of motion are written in terms of modal co-ordinates. However, when the system components experience non-linear deformations the use of reduction methods is not possible, in general, and the ÿnite element nodal co-ordinates are the generalized co-ordinates used. Furthermore, depending on the type of ÿnite elements used to represent each exible body, the nodal co-ordinates may include all node rotations and translations or only some of each. Regardless of the type of generalized co-ordinates adopted it is required that kinematic joints are deÿned. The complete set of joints available in a general-purpose multibody code must include, for each particular type of joint, restrictions involving only rigid bodies, or only exible bodies, or exible and rigid bodies. Therefore, the e ort put into the development and implementation of any joint is at least three times as much as the initial work done in the implementation of joints with rigid bodies only. The concept of virtual bodies provides a general framework to develop general kinematic joints for exible multibody systems with minimal e ort, regardless of the exible co-ordinates used. Initially, only a rigid constraint between the exible and a massless rigid body is developed. Then, any kinematic joint that involves a exible body is set with the massless rigid body instead, using the regular joint library of the multibody code. The major drawback is that for each kinematic joint involving a exible body it is required to use six more co-ordinates per virtual body and six more kinematic constraints. It is shown in this work that for small elastic deformations, for which the mode component synthesis is applied, the use of sparse matrix solvers can compensate for the computational overhead of involving more co-ordinates and kinematic constraints in the system, due to the use of virtual bodies. For non-linear deformations, where the generalized co-ordinates are the global positions of the ÿnite-element nodes, the use of the virtual body concept does not require an increase in the number of system co-ordinates or kinematic constraints. By introducing the rigid joint between the exible body nodal co-ordinates and the virtual body, with the use of Lagrange multipliers, and then solving the equations explicitly for * Correspondence to: Jorge Ambrà osio, Instituto de Engenharia Mecânica (IDMEC), Instituto Superior Tà ecnico, Av.Rovisco Pais, 1049-001 Lisboa, Portugal. † E-mail: jorge@lemac.ist.utl.pt these multipliers the resulting equations of motion for the subsystem have the same degrees of freedom as the original exible body alone. The di erence is that degrees of freedom associated to the virtual body are used as co-ordi...