The simulation of flexible multibody systems with unilateral contact conditions and impacts requires advanced numerical methods. The nonsmooth generalized-α method was developed in order to combine an accurate and second-order time discretization of the smoother part of the dynamics and a consistent but first-order time discretization of the impulsive contributions. Compared to the Moreau-Jean scheme, this approach improves the quality of the numerical solution especially for the representation of the vibrating response of flexible bodies. It entirely relies on the formal definition of a so-called smooth motion that captures a non-impulsive part of the total nonsmooth motion. This definition may account for some contributions of the bilateral constraints and/or of the active unilateral constraints at velocity or at acceleration level. This chapter shows that the formulation of the constraints strongly influences the numerical stability and the computational cost of the method. A strategy to enforce the bilateral and unilateral constraints simultaneously at position, velocity and acceleration levels is also established with a careful formulation of the activation criteria based on augmented Lagrange mul-tipliers. In the special case of smooth systems, a comparison is made with more standard solvers for differential-algebraic equations. The properties of this method are demonstrated using illustrative numerical examples of smooth and nonsmooth mechanical systems.