Flexible satellite structures (e.g., solar arrays, booms) demonstrate complex nonlinear dynamics arising from structural/material nonlinearities, environmental disturbances, and control interactions. Key phenomena include internal resonance—occurring when natural frequencies become commensurate, enabling strong intermodal energy transfer—and limit cycle oscillations (LCOs), self-sustained vibrations induced by geometric/thermal nonlinearities. These effects are particularly pronounced during eclipse transitions, where transient thermal loads interact with vibration modes, potentially amplifying responses. Lightweight multi-DOF systems are especially vulnerable, as internal resonance can exacerbate vibration amplitudes while LCOs persist without external excitation, collectively challenging stability and performance in space missions.

Thermally induced vibrations in satellite structures have been extensively studied through various analytical and numerical approaches. Gulick and Thornton [1] pioneered the analysis of thermal-induced motion in spin-stabilized satellites, employing Floquet theory to identify stability boundaries. Subsequent studies have expanded this foundation through coupled thermo-structural analyses. Duan et al. [2] developed geometric nonlinear methods for large space frames, while Shen and Hu [3] characterized thermal flutter mechanisms in cantilevered solar panels. Zhang et al. [4] established stability criteria for thermally excited beams, providing fundamental insights into nonlinear dynamic responses.Recent advances in modeling techniques have enabled more sophisticated investigations. Shen et al. [5] formulated an absolute nodal coordinate-based unit for analyzing large flexible structures, demonstrating its utility as a modular component for complex systems. Their subsequent work [6] quantified thermal shock effects on deploying booms, revealing transient response characteristics. Wei et al. [7] contributed a global mode method (GMM) for flexible manipulators, enhancing dynamic modeling accuracy for multi-body systems.Comparative studies have validated analytical predictions against numerical benchmarks. Parisse and Angeletti [8] developed thermal shock vibration models for slender spacecraft with solar panels, verifying results against finite element analyses. Lu et al. [9] advanced coupled thermal-structural-attitude dynamics for flexible beams, demonstrating how orbital motion influences flutter stability under thermal loading. Zhang et al.[10] extended this to solar sails, identifying coupled radiation-thermal-structural interactions that govern flutter behavior.Specialized analyses have addressed distinct structural configurations. Shen and Hu [11] employed thermoelastic-structural frameworks to analyze thin-walled beams under solar flux shocks. Qian et al. [12] discovered novel coupled-direction flutter modes in solar sails using traveling wave analysis. Li and Yan [13] optimized composite honeycomb arrays by characterizing vibration responses to orbital thermal variations.Nonlinear phenomena have received particular attention. Li et al. [14] established finite element-based flutter criteria for thin-walled structures during orbital transitions. Yuan and Xiang [15] derived Lyapunov-stable solutions for cantilevered beams under solar heating, validating theoretical predictions numerically.Application-focused studies have provided critical design insights. Thornton and Kim [16] developed coupled models for rolled-up solar arrays, identifying key dimensionless parameters through Hubble Telescope case studies. Thornton et al [17] further analyzed self-shadowed arrays, demonstrating how parallel-member configurations induce stable vibrations.Rigid-flexible coupling effects have been systematically investigated. Yuteng Cao, et al [18–20] developed hierarchical modeling approaches, employing Hamiltonian mechanics and GMMs to capture maneuvering-induced vibrations.He et al [21] created efficient low-dimensional models for array-hinged spacecraft using Rayleigh-Ritz approximations. Liu et al. [22] revealed attitude control challenges from thermally excited vibrations.Broader dynamic phenomena have also been explored. Aslanov [23] characterized chaotic attitude dynamics in LEO satellites with flexible panels, identifying critical altitudes for aerodynamic stabilization. Gao et al. [24] analyzed internal resonance effects in space antennas, deriving nonlinear frequency-amplitude relationships through multiscale methods.Collectively, these studies establish comprehensive frameworks for understanding thermally induced vibrations across satellite subsystems, though gaps remain in predicting coupled thermo-structural-control instabilities during complex orbital maneuvers. Motaharifard et al.[25] a nonlinear analytical model of a spacecraft with flexible solar arrays, investigating geometric nonlinearities and internal resonance phenomena during the penumbra phase using the method of multiple scales (MMS). Ji and Zhang[26] conducted a dynamic analysis of rigid–flexible coupling spacecraft based on Euler parameters

The Hubble Space Telescope's flexible solar arrays exhibit complex nonlinear dynamics when subjected to thermal gradients during penumbra transitions. This study investigates thermally induced LCOs and internal resonance phenomena arising from rigid–flexible–thermal coupling. Using the method of multiple scales, we identify critical internal resonance conditions governing modal interactions. The analysis employs a hyperbolic heat conduction model with second-order time derivatives to capture finite thermal wave propagation speeds (particularly at 0.95c) and material response delays during rapid thermal transients. Nonlinear cubic terms model radiation flux separation during eclipse transitions, while the Wagner function characterizes transient thermoelastic responses.Combined analytical-numerical results reveal how structural nonlinearities (from large deformations) and thermal nonlinearities (from time-varying heat loads) collectively drive LCO formation. Phase portraits, Poincaré maps, and bifurcation diagrams demonstrate period-doubling routes to instability in the thermal flutter regime. Key findings include: (1) the critical role of thermal wave velocity ratios near 0.95c in triggering coupled-mode oscillations, (2) Wagner-function-predicted transient amplitudes matching numerical simulations within 7%, and (3) identification of three distinct LCO regimes governed by the interplay between structural stiffness and thermal relaxation times. These results provide new insights for mitigating performance-degrading vibrations in precision space telescopes undergoing orbital thermal transitions.

This graph illustrates the changes in the speed ratio of thermal wave propagation relative to the speed of light as a geostationary satellite experiences varying conditions of shadow and penumbra during a solar eclipse. The time interval from 90 to 120 seconds shows a gradual decrease in the speed ratio from 1 to 0, indicating the satellite's transition into the shadow [33]. From 120 to 150 seconds, the speed ratio remains consistently at 0, reflecting that the satellite is fully engulfed in the shadow. Subsequently, between 150 and 180 seconds, the speed ratio gradually increases from 0 to 1, marking the satellite's exit from the shadow. Following this, from 180 to 240 seconds, the satellite returns to normal conditions, with the speed ratio stabilizing at 1. The distinct regions of shadow, penumbra, and the respective entry and exit times are clearly highlighted in the graph using different colors, facilitating an easy understanding of the variations in speed ratio throughout the eclipse phases [34].

To investigate internal resonance and LCOs in a flexible solar array model, representative physical parameters have been adopted based on previous studies [19]. The thermal properties responsible for inducing structural vibrations are summarized in Table 1, as reported in [35].

For validation purposes, the results derived from the FEM were benchmarked against those reported in the reference study, demonstrating excellent consistency. Furthermore, the GMM was utilized to perform modal analysis and extract mode shapes. This approach proves highly effective in capturing the dynamic characteristics of the system under combined thermal and mechanical loading conditions. The strong correlation between the outcomes of both methods and those of the reference confirms the accuracy and robustness of the proposed model. Table 3 and Fig. 4 illustrate this comparison. In this context, the variables u, w, and α denote the in-plane bending, out-of-plane bending, and torsional responses, respectively.

In the coupled structural-thermal analysis of satellites, if the in-plane bending stiffness is sufficiently high, its effects can be neglected. This assumption simplifies the nonlinear BBT equations to the more manageable bending-torsion (BT) equations, which focus on the interaction between these two dominant deformation modes. By adopting this simplification, the structural-thermal analysis becomes more efficient, allowing for a clearer investigation of the primary effects of bending and torsion. The presented equations are formulated in a dimensionless form, ensuring a direct correlation between the structural properties and the thermal response of the system.

In a linear thermoelastic system, the dynamic characteristics can be distinguished based on the relative speed of thermal wave propagation compared to the flutter velocity (VF). In other words, the rate at which heat propagates through the structure can significantly influence its stability and dynamic behavior. If the thermal wave propagation speed exceeds the flutter velocity, temperature variations rapidly spread throughout the structure, potentially inducing swift changes in its mechanical properties and compromising its stability. Conversely, if the thermal wave propagation speed is lower than the flutter velocity, temperature changes disperse more gradually, reducing their impact on the system’s dynamics. These distinctions play a crucial role in the precise analysis of thermal flutter in linear thermoelastic systems and can significantly affect the likelihood of flutter occurrence[45].

Table 1, 2 presents the thermal and structrual characteristics [45]. In the case where the ratio of the thermal wave propagation speed to the speed of light is approximately 0.95, the system damping reaches zero, as shown in Fig. 5. This condition is referred to as the thermal flutter speed range. When the velocity is less than VF​ (Vr​=0.85), the responses decay, as shown in Fig. 6. The possibility of a 3:1 internal resonance and the associated response characteristics are investigated. As shown in Fig. 5, the maximum frequency ratio ω/ωα is 4.6, with the peak occurring at V = 0. Since this ratio exceeds 3 and ω > ωα internal resonance does not occur in the structural dynamics of the flexible arrays.

Eventually, the responses decay to zero. For velocities greater than VF​, the responses increase over time, as illustrated in Fig. 7.

As illustrated in Fig. 5, the coupled thermal-structural frequencies corresponding to the bending and torsional modes differ under the velocity conditions Vr​=0.85 and Vr​=1.05. However, the responses indicate synchronization between the two modes, which represents a characteristic feature of coupled thermal-structural behavior. Due to the presence of two distinct types of nonlinearities and their possible interactions, three representative cases are examined, each corresponding to a different combination of these nonlinear effects. To investigate the existence of LCO at velocities below VF​, a reference case is considered at the velocity condition Vr​=0.95. Since closed orbits in phase portraits indicate periodic motions, phase portraits are employed to verify the existence of LCOs. The response characteristics are summarized in Table I. In this table, the letters ‘W’ and ‘T’ are used to denote the types of nonlinearities considered in each case. The case labeled ‘BT3’, which includes both ‘W’ and ‘T’, indicates that the system exhibits structural nonlinearities due to solar array flexibility (‘W’) as well as thermal nonlinearity (‘T’). As shown in the table LCO type responses are observed in cases that involve structural nonlinearities.

For velocities below the flutter threshold, the motion decays over time, as illustrated in Fig. 6, whereas for velocities above the flutter point, the response exhibits growth, as shown in Fig. 7. At Vr​=0.95, the response characteristics of certain cases involving structural nonlinearities exhibit strong similarities. Therefore, phase portraits and time histories for a representative case are presented in Fig. 8. The time histories demonstrate that the response amplitudes remain sustained over time. The corresponding phase portraits and Poincaré maps are presented in Fig. 9. In this figure, the phase portrait is represented by a solid line, while the Poincaré maps are indicated by the symbol ‘O’. The presence of an isolated closed orbit in the phase portrait confirms the existence of a limit cycle. From the phase portrait, an isolated closed orbit can be identified. This closed orbit corresponds to a periodic motion[43], as illustrated in Fig. 8. This type of periodic motion is characterized by two fixed points, AF and BF. These figures indicate that for any case involving structural nonlinearities, LCO type responses can occur at velocities below VF​. Although phase portraits offer valuable insights into LCO-type responses, they are not well-suited for summarizing key characteristics. To examine the effects of velocity variation and the relative motion between bending deflection and torsional response, Poincaré maps and bifurcation diagrams are employed. For a forced vibration system, a Poincaré map can be constructed using a sampling rule that aligns with the period of external excitation. This approach allows for the analysis of the system’s long-term behavior by capturing its state at discrete, periodic intervals.

Since steady-state responses can occur in thermoelastic systems even in the absence of external excitation, alternative sampling rules are necessary for constructing Poincaré maps in such cases. To construct the bifurcation diagram of pitch amplitude and the Poincaré map of pitch motion, a systematic sampling of the system’s dynamic response is performed over a range of velocity conditions, allowing the identification of qualitative changes and the onset of nonlinear phenomena such as LCO or bifurcations [46, 47]. Two sampling rules are employed: the condition ​=0, indicating zero velocity of the pitch angle, and h = 0.0, corresponding to zero bending deflection, respectively. These rules are used to construct the Poincaré map and capture the system’s dynamic behavior at specific, meaningful states. In this study, an additional sampling rule, , is adopted to generate Poincaré maps of the pitch motion. This condition enables the capture of the system’s dynamic behavior at instances when the rate of change of the bending displacement is zero, providing further insight into the nonlinear response characteristics. By applying this new sampling rule, equilibrium points where both and can be identified in the Poincaré map. If the fixed points in the Poincaré map satisfy the condition , these points correspond to equilibrium states of the system. By utilizing the fixed points observed in the Poincaré maps for various velocity conditions, bifurcation diagrams are generated.

These diagrams provide a clear visualization of the system's dynamic transitions, revealing the relationship between system behavior and varying velocity parameters. The bifurcation diagram for a case with only thermal nonlinearities is shown in Fig. 10. This diagram illustrates the system's dynamic behavior and the critical transitions associated with thermal effects under varying conditions. As illustrated in Fig. 10, the number of fixed points varies with the velocity conditions: a single fixed point is observed Vr​<1, while two fixed points appear when Vr​>1. This indicates a bifurcation behavior influenced by the system's velocity. As illustrated in the phase plane, the motion in the upper half-plane (where ) must progress to the right, since the variable on the horizontal axis (ϕ) increases with time when . Conversely, for , the motion proceeds to the left. This interpretation holds under the assumption that ​. However, if the convention ​ is adopted, the direction of motion must be reversed accordingly.

This indicates a bifurcation behavior influenced by the system's velocity. The single fixed point in the region of Vr​<1 indicates that the responses decay to zero, as shown in Fig. 6. However, the two fixed points at velocities greater than VF ​suggest the presence of LCOs, as demonstrated in Figs. 11 and 12. As shown in Fig. 12, the corresponding phase portraits and Poincaré map reveal closed orbit(s) and fixed points, respectively, indicating periodic motion in the system.. For a velocity condition of Vr​<1, the responses decay to zero in a manner similar to that shown in Fig. 6. In contrast, for a velocity condition of Vr​>1, the responses grow in a manner similar to Fig. 7.

This indicates a transition from a stable response to an unstable one as the flexible satellite structure moves out of the shadow. Based on the results presented in Table 4, the effects of velocity variation are investigated for selected cases involving structural nonlinearities. Initially, a case featuring only structural nonlinearity is examined to isolate and understand its influence on the system's dynamic response.

The bifurcation diagram for a case in which both forms of nonlinearity—structural and thermal—are considered is presented in Fig. 13. This figure illustrates the complex dynamic behavior resulting from the interaction between these nonlinear effects. As shown in this figure, multiple fixed points with varying magnitudes of the pitch angle amplitude (α) are observed for Vr​<1. This indicates the presence of multiple steady-state solutions and highlights the complexity of the system's nonlinear dynamic behavior under subcritical velocity conditions. These fixed points indicate that LCO type responses are possible within the subcritical velocity range Vr​<1. This suggests the existence of sustained periodic motion even before reaching the flutter velocity, due to the nonlinear characteristics of the system. For Vr​<VC​, the responses decay to zero.

In the intermediate velocity range VC < Vr​<VF​, at least three fixed points exist, including one at zero pitch amplitude (α = 0). The presence of this zero fixed point suggests that the system's response may decay, depending on the initial conditions, indicating potential coexistence of stable and unstable solutions within this velocity range.

For small initial conditions, the response tends to decay to zero; however, larger initial conditions may lead the system to converge toward other fixed points with non-zero amplitudes. This behavior reflects the sensitivity of the system's dynamic response to initial states and indicates the presence of multiple attractors in the nonlinear regime. This type of subcritical bifurcation behavior persists within this velocity range, indicating the coexistence of multiple steady-state solutions and the potential for abrupt transitions in response depending on initial conditions.

The bifurcation diagram can be partitioned into several distinct regions, each representing different qualitative behaviors of the system's dynamic response as a function of velocity. Three fixed points are observed in regions B and D. In region C, the number of fixed points exceeds three. Although the number of fixed points in regions B and D is the same, their corresponding amplitudes differ. This highlights the complexity of the system's behavior within these regions. As illustrated in Figs. 14 and 15, these regions exhibit distinct dynamic characteristics.

The Poincaré maps provided in these figures clearly reveal the positions of the fixed points labeled AF​, BF​, and CF​, facilitating the identification and analysis of the system's qualitative behavior in each region. The two fixed points, AF​ and BF​, correspond to periodic motions. The presence of a single AF​ (or BF​) indicates that the bending and torsional motions share the same (nonlinear) period, reflecting synchronized dynamic behavior between the two modes. The third fixed point, CF​, corresponds to a damped response. Depending on the initial conditions, the system's behavior either decays to zero or exhibits LCO, indicating the presence of a subcritical bifurcation phenomenon.

The two fixed points, AF​ and BF​, presented in Fig. 15(a), are located near the condition ; however, the remaining fixed points in the same figure are situated closer to the maximum value of ​, highlighting distinct dynamic behaviors within the system. Because the sampling condition is , the fixed points located near zero indicate that the relative phase angle between the bending and torsional motions is also close to zero, as illustrated in Fig. 15(b). However, in the other case where is near its maximum, the two motions exhibit an approximate 90-degree phase difference, as illustrated in Fig. 15(b).

Although the bending motion is not explicitly presented in these figures, its relative behavior can be inferred from the Poincaré map. In region C, two intriguing types of responses are observed, as illustrated in Figs. 16, 17, and 18. At a velocity near Vd​, the nonlinear frequency ratio approaches 3, as illustrated in Fig. 16. In this figure, three points corresponding to AF​ and three points corresponding to BF​ are identified. This indicates that the period of the bending motion is approximately three times that of the torsional motion. For another velocity close to Ve​, quasi-periodic motion is observed, as illustrated in Fig. 18. As illustrated in Fig. 18, the fixed points form a circular pattern. This type of Poincaré map typically emerges when the frequency ratio between bending and torsional motions is irrational, indicating quasi-periodic behavior in the system.

A nonlinear structural model of a flexible satellite is developed based on the extended Hamilton's principle. The BBT equations of motion are formulated up to third-order nonlinear terms. The equations incorporate the nonlinearity induced by heat flux variations resulting from the separation of radiant energy from the solar array surfaces during penumbra, to accurately capture the thermal-structural coupling effects. In a scenario where the satellite operates in space and is excited exclusively by solar radiation, the structural frequencies are influenced by the propagation speed of the incident solar waves. Under such conditions, the excitation frequencies of the structure can exhibit a strong dependence on the wave propagation velocity, potentially altering the dynamic response of the system. Internal resonance conditions are identified using the Method of Multiple Scales. Two resonance cases, ​and , arise due to thermal nonlinearities, while two additional conditions ​, and ​ result from structural nonlinearities. Among these, two resonance conditions and are associated with cubic nonlinearities, while the other two ​ and arise from quadratic nonlinearities. These internal resonance phenomena are investigated in the context of a solar array mounted on a flexible satellite. The conditions for potential internal resonance vary based on factors like the stiffness ratio, , the length of the solar arrays, and thermal effects. The solar array model considered in this study does not satisfy the resonance condition, as the maximum frequency ratio occurs at V = 0.Possible causes of LCO responses in flexible satellites are investigated using solar array models inspired by the flexible arrays of the Hubble Space Telescope (HST). For a specific case, the effect of each individual nonlinearity is examined in detail. Bifurcation diagrams are used to demonstrate the effect of changes in velocity. The bifurcation diagrams are constructed using Poincaré maps. For the system exhibiting thermal nonlinearity, the bifurcation diagram reveals a supercritical bifurcation behavior. Stable solutions are observed in both velocity regions, Vr < 1 and Vr > 1. In the case of structural nonlinearity in the flexible arrays, a subcritical bifurcation is identified. Depending on the initial conditions, the system's behavior either decays to zero or exhibits LCO, indicating the presence of a subcritical bifurcation phenomenon. The bifurcation diagram can be partitioned into several distinct regions, each representing different qualitative behaviors of the system's dynamic response as a function of velocity. For small initial conditions, the response tends to decay to zero; however, larger initial conditions may lead the system to converge toward other fixed points with non-zero amplitudes. This behavior reflects the sensitivity of the system's dynamic response to initial states and indicates the presence of multiple attractors in the nonlinear regime. Results indicate that structural nonlinearities significantly amplify oscillation amplitudes, while thermal nonlinearities contribute to asymmetric response patterns and increased sensitivity to initial conditions. Parametric studies show that even slight changes in system properties such as mass or stiffness can shift critical velocities and fundamentally alter the dynamic response. This research underscores the critical role of nonlinear thermal-structural interactions in the dynamic behavior of flexible space systems and offers a comprehensive analytical framework to guide the design and safe operation of future spacecraft exposed to severe thermal and mechanical environments.It is also important to note For model validation, the results obtained through the FEM were compared with those from the reference paper, showing excellent agreement. Additionally, the GMM was employed for mode analysis and to extract the modal shapes.Implementing control strategies such as active damping mechanisms, adaptive structural tuning, and smart materials can help minimize unwanted vibrations. Techniques like vibration absorbers, real-time feedback control, and optimized structural configurations can also enhance stability. Future research should focus on integrating these strategies into spacecraft design to ensure reliability and performance in space environments.