Finite element methods are widely used for simulations of contact dynamics of flexible multibody systems. Such a simulation is computationally very inefficient because the system’s dimension is usually very large and the simulation time step has to be very small in order to ensure numerical stability. A potential solution to the problem is to apply a model reduction method in the simulation. Although many model reduction techniques have been developed, most of them cannot be readily applied due to the high nonlinearity of the involved contact dynamics model. This paper presents a solution to the problem. The approach is based on a modified Lyapunov balanced truncation method. A numerical example is presented to demonstrate that, by applying the proposed model reduction method, the simulation process can be significantly speeded up while the resulting error caused by the model reduction is still within an acceptable level.