© 2016 Elsevier Ltd. All rights reserved. A computationally efficient method is presented for static and dynamic analysis of a beam-like structure on a viscoelastic foundation with the unilateral contact constraint at their interface. The non-smooth dynamics of the system is modelled using the Euler-Bernoulli theory for the beam and the bilinear Winkler model for the substrate. Thus, the unilateral contact is the only source of non-linearity. The proposed approach relies on reducing the non-smooth continuous system to a piecewise smooth multi-degree of freedom model. The beam is represented as a chain of discrete units through the use of the lattice spring model (LSM) and consequently, the phase space is divided into a number of subdomains. Hence, the system smoothness is lost when the trajectory of lattice nodes cross the boundaries between these subdomains. An effective algorithm has been developed to handle the unilateral constraints by tracking the trajectories and successively capturing the intervals that the nodes spend in each of the smooth regions. Unlike more commonly used methods, it neither relies on any prior information (such as the number and location of the lift-off areas) nor on the non-linear solvers. The accuracy and robustness of this numerical technique have been demonstrated in several application examples. Furthermore, the developed algorithm is utilised in combination with the shooting method and continuation of the periodic motions to obtain the non-linear normal modes (NNMs) of the system. We show that frequency-stiffness (Ωi u-k̄w) plots can be used to describe the salient features of the NNMs. For the beam on the full foundation, the first mode in the configuration space has a linear shape and thus, the conventional bilinear formulation can accurately approximate the non-linear frequency. For higher modes, the modal lines are open curves and the bilinear approximation loses its accuracy. In addition, internal resonances may occur once the frequencies are commensurate. They are characterised as new multimodal motions emanating from the backbone in the frequency-stiffness plot.