The increasing reliance of aerospace structures on numerical analyses encourages the development of accurate, yet computationally efficient, models. Finite element (FE) beam models have, in particular, become widely-used approximations during preliminary design stages and to investigate novel concepts, e.g. aeroelastic tailoring. Over the last 50 years, developments in hp-FE methods based on elements of variable size (h) and polynomial degree (p) have helped reduce the computational cost of numerical analyses. Concurrently, many structures, including aircraft wings and wind turbine blades, have gradually increased in length, slenderness, and complexity. As a result, modern blades and wings regularly operate beyond the range of linear deformation, requiring non-linear analyses for which hp-FE methods are often not readily applicable. The aim of this paper is, therefore, to derive a co-rotational finite element formulation for enriched 3-, 4-and 5-noded beam elements, suitable for non-linear hp-FE refinement. To this end, we derive the mathematical formulation to incorporate enriched elements within the co-rotational FE beam framework. The proposed formulation is then validated against multiple non-linear benchmark problems and an experimental case study.