Flexible beams are becoming ubiquitous in several industrial applications, as new projects often aim for lighter and longer structures. This fact is directly related to the new challenging demands on structural performances, or it is a simple consequence of the engagement of industries in cost reduction programs (usage of less material). Flexible beams are usually modeled under the assumption of large displacements, finite rotations, but with small strains. Such hypotheses allow the equation of motion to be built using co-rotational finite elements. The co-rotational formulation decomposes the total motion of a flexible structure into two parts: a rigid body displacement and an elastic (small) deformation. This way, the geometric nonlinearity caused by the large displacements and rotations of the beam s cross sections can be efficiently computed. One of the novelties of this thesis is the direct usage of the equation of motion generated by a co-rotational finite element formulation in the computation of nonlinear normal modes (NNM). So far, most of the dynamic analyses with co-rotation finite element models were restricted to numerical integrations of the equation of motion. The knowledge of NNMs can be beneficial in the analysis of any nonlinear structure since it allows a thoroughly understanding of the vibratory response in the nonlinear regime. They can be used, for example, to predict a hardening/softening behavior, a localization of the responses, the interactions between modes, the existence of isolas, etc. The Rosenberg s definition of NNM as periodic solutions (non-necessarily synchronous motion) is adopted here. The Harmonic Balance method and the Shooting methods are presented and used to compute periodic solutions of nonlinear systems. A numerical path continuation scheme is implemented to efficiently compute NNMs at different energy levels. Numerical examples show the capability of the proposed method when applied to co-rotational beam elements.