Nonlinear phenomena are becoming increasingly significant in modern mechanical systems due to the transition to lighter, slender structures, and the will to exploit nonlinearities to design more efficient systems. This shift requires accurate system identification techniques to develop reliable models, which are crucial for both the design and monitoring of these systems. However, the identification process is computationally intensive, as numerical models can reach millions of degrees of freedom, with nonlinearities adding significant complexity to the computations. This work aims to address a gap in the current literature related to the simulation of periodic structures with nonlinear boundaries. We propose an efficient computational framework that not only reduces the size of the problem but also provides more physical insights than a classical finite element approach. The proposed method is an extension of the Wave Finite Element Method (WFEM), which combines the Floquet-Bloch theory and the Finite Element Method, to compute the dynamics of periodic waveguides with nonlinear boundaries. The Harmonic Balance Method is used to recast the governing equations into a nonlinear algebraic system, which can then be solved using numerical continuation algorithms. To demonstrate the validity and benefit of this approach, the predictions derived from the proposed methodology are compared to results obtained through conventional Finite Element analysis. Excellent agreement and a notable reduction in computational cost are observed. This methodology enables simulations of realistic models within a feasible timeframe. This aligns with the demands of an R&D product cycle, which requires simulation intensive tools for system identification and digital twinning.