Fast and accurate simulation of mechanical structures with complex geometry requires application of the finite element method. Τhis leads frequently to models with a relatively large number of degrees of freedom, which may also possess nonlinear properties. Things become more complicated for systems involving unilateral contact and friction. In classical structural dynamics approaches, such constraints are usually modeled by special contact elements. The characteristics of these elements must be selected in a delicate way, but even so the success of these methods can not be guaranteed. This study presents a numerical methodology, which is suitable for determining dynamic response of large scale finite element models of mechanical systems with multiple unilateral constraints. The method developed is based on a proper combination of results from two classes of direct integration methodologies. The first one includes standard methods employed in determining dynamic response of structural models possessing smooth nonlinearities. The second class of methods includes specialized methodologies that simulate response of dynamical systems with unilateral constraints. The validity and effectiveness of the methodology developed is illustrated by numerical results.