## 1 Introduction  
In recent decades, the increased requirements of mechanical and physical performances in sustainable structures have favored the use of carbon nanotubes (CNTs) as reinforcement phase in composite materials, combined with the concept of functional graduation of materials to yield smooth and continuous characteristics [1, 2, 3, 4, 5]. Among different composite reinforcements, CNTs offer outstanding mechanical properties in terms of stiffness, strength, and mass density. The addition of only 1% of CNTs in the composite material,for example, provides an increased elasticity modulus of about 40%, with a ductility reduction, an increased fracture toughness [6], and improved damping properties [7]. Among composite applications, like rotocrafts, robotic manipulators, or wind turbine blades, many structural components, like beams, plates, and shells, can be subjected to aerodynamic loadings, with the introduction of structural instabilities, e.g., flutter and divergence, with possible catastrophic effects. In most cases, the vibration amplitudes are large, making various forms of nonlinearity significant and situating the problem within the domain of nonlinear structural dynamics [8]. The dynamic behavior of flapping beams, for example, is expected to exhibit many complex phenomena commonly associated with nonlinear systems--such as bifurcations, quasi-periodicity, and chaos [8; 9; 10] and more specifically, those seen in nonlinear beam dynamics [11; 12; 13; 14; 15; 16]. In such a context, in recent decades, an increased attention has been devoted to understanding the nonlinear structural dynamics through both numerical simulations and experimental studies. The mathematical and computational models employed in these investigations span a wide range of complexity--from finite element models of geometrically exact, curved composite beams [17] to simplified, single-mode representations capturing chaotic vibrations in beams with nonlinear boundary conditions [13]. The research on nonlinear dynamics, and particularly in the area of beam dynamics, is quite extensive, but more comprehensive discussions can be found in references [11; 18] for general nonlinear dynamics, [9] for nonlinear vibrations, and [19] for nonlinear beam vibrations. It is well known from the literature that nonlinear systems can exhibit a wide variety of resonance phenomena, influenced by factors such as the strength of nonlinearity, the system's degrees of freedom, and the characteristics of the external excitation [20]. Thus, extensive studies have been conducted on the primary and superharmonic resonance behaviors of fundamental nonlinear systems, such as the Duffing oscillator, with particular attention to stability characteristics and the influence of key system parameters [21; 22; 23]. [PERSON] and [PERSON] [22], for instance, compared second-order multiple scales perturbation techniques with numerical time integration to analyze the steady-state response and stability near the third-order superharmonic resonance. Their findings indicated that second-order perturbation results closely mirrored those from the first-order analysis and only modestly aligned with numerical outcomes. More recently, [PERSON] et al. [24] employed a time-spectral approach to investigate the complex response behavior of the Duffing system, particularly focusing on the third-order superharmonic regime. Their analysis detailed the emergence of multi-valued response branches and jump phenomena, examined how harmonic components varied with frequency, and evaluated the roles of damping and excitation amplitude on the overall response landscape. In addition to these theoretical investigations, different computational strategies have been increasingly employed in the literature to study the aeroelastic behavior of various structural members. Among the recent works, a classical finite approach was proposed in refs.