1. Introduction
While orbiting the Earth, satellites undergo periodic heating and cooling due to alternating exposure to sunlight and shadow, resulting in varying thermal conditions. Rapid changes in temperature gradients can induce relevant vibrations in the spacecraft’s flexible appendages (e.g., solar panels and antennas) due to their low stiffness, mass and heat capacity. At the same time, the elastic deformation of in-orbit structures is coupled with the rigid-body dynamics of the main hosting platform, thus leading to disturbances in the spacecraft’s attitude [
1,
2,
3]. Indeed, thermally induced vibrations are generally observed indirectly by analyzing telemetry data from on-board attitude sensors, and they are rarely detected using a dedicated measuring setup, as in the case of in situ reflectors on solar panels [
4], due to related additional mass and costs. Notable research efforts are, hence, addressed at studying and predicting such a phenomenon, as it can affect stability and pointing performance or even lead to structural damage. Historic examples of thermo-elastic–rigid coupling in spacecraft dynamics include the Alouette I and Explorer XX satellites, which experienced anomalies in their spin behavior due the joint action of solar heating-induced bending and radiation pressure in long tubular beams [
5]. An attitude deviation from the local vertical was observed in gravity gradient-stabilized satellites equipped with long booms, such as the 1963-22A [
6]. This was traced back to vibrations of the 30-m long, slender structure realized in beryllium–copper material, which was coated with silver in the subsequent satellites to successfully reduce the phenomenon. Impulsive disturbance torques were noticed during the operations of the TOPEX satellite and then assessed as the consequence of the fast heating of its solar array during the sun–shadow transition [
7]. Later, a famous case study involved the investigation of orbital transition disturbances in the Hubble Space Telescope, which were then associated with the eclipse–sunlight transition and deformation of its solar array panel [
8].
Thermally induced vibrations were theoretically predicted before being observed in orbit by Boley [
9], who theorized that, when a beam undergoes rapid external heating, a time-dependent thermal moment develops due to the temperature gradient between its heated and unheated sides. By introducing a basic adimensional parameter, Boley demonstrated that, if the thermal response time is comparable to or shorter than the structural response time, transient vibrations occur, which is generally the case for lightweight and flexible space structures. The study was then extended to more complex structures, such as plates [
10], and an expression of the ratio between the maximum dynamic and static deflection was found [
11]. Later, an analytical formulation to simulate the thermo-structural behavior of space flexible panels was derived by Thornton et al. [
12], also identifying a criterion to detect thermal flutter. Moreover, Frisch [
13], Tsuchiya [
14] and Li et al. [
15] addressed the problem of the influence of thermally induced dynamics on the attitude (i.e., orientation) of satellites. Numerical and analytical methods were verified via experimental studies, such as the ones conducted by Johnston and Thornton [
16] and Su et al. [
17] to reproduce thermal snap and induced vibrations of space structures, respectively.
As space structures have grown larger and more complex, numerous researchers have focused on developing advanced numerical and analytical methods to study thermally induced vibrations. In this scenario, the finite element method (FEM) was first proposed by Mason [
18], and it is now a standard procedure for advanced the design stages of spacecraft systems. However, simulations are generally time-consuming due to the thermal–structural coupling effect and dynamic nonlinearities [
19], and usually, FEM models are developed for structural dynamics simulation and are not tailored to thermo-mechanics analyses [
20]. Solutions have been developed to mitigate this issue, such as the application of Fourier FEM by Xue et al. [
21], which was only recently experimentally tested in a laboratory environment by Wang et al. [
22]. A different approach using absolute nodal coordinate formulation was proposed by Shen et al. [
23] for composite solar panels and axially moving deployable structures [
24]. Mu et al. [
25] developed in-orbit governing equations via Hamilton’s principle for a large, flexible beam structure under the influence of a gravity gradient, solar radiation pressure and thermal radiation, discussing system instability due to the resonance between the attitude motion and the vibration and identifying the thermal input as the dominant effect on structural vibrations in GEO orbits. A methodoloy based on the Global Mode Method (GMM) and Finite Difference Method (FDM) was developed by Cao et al. [
26] to simulate thermal alternation-induced vibrations for solar arrays. Thermal–structural analyses for a flexible spacecraft with double solar panels were established by Liu et al. [
27] via a two-way coupled Hamiltonian principle and finite difference formulation to solve the transient heat-conduction problem. An analytical-numerical solution to the thermo-elastic problem in aluminum-composite panels during an eclipse was developed by Ganilova et al. [
28] to reproduce the experimental results obtained in [
29].
Over the years, several strategies have been proposed to address thermally induced vibrations and model coupling effects on spacecraft under heating and cooling conditions. However, these approaches often result in complex and computationally expensive models that, while elegant, are not easily accessible for quick estimations or a preliminary engineering analysis of the phenomenon driving the system dynamics. This paper addresses the need for a more straightforward yet effective analytical method to capture the essential physics of the coupled rigid–flexible–thermal dynamics of spacecraft, while being meaningful for engineering purposes. Specifically, it focuses on long, slender components, which are particularly susceptible to disturbances due to differential thermal stresses. The core novelty of this study lies in the proposed modeling method to describe the coupling between thermal and structural responses and capture the potential thermal flutter behavior of orbiting beam-like structures. Drawing on the fact that the first bending elastic mode is the most excited due to differential temperatures in the two opposed surfaces, the paper proposes a method that conceptualizes the beam-like structure as two halves connected via a cylindrical hinge, reproducing a deformation similar to the first modal shape. This abstraction transforms the thermal–structural coupling problem into a rotational mass–spring–damper system. The parameters of this system—such as rotational stiffness, damping and mass—are derived from physical considerations, making the model both practical and physically grounded. The results indicate that such an approach provides sufficient accuracy when compared to an equivalent finite element model in predicting thermal jitter issues while significantly reducing computational demands and modeling efforts for spacecraft early design phases. The significance of this research lies in its ability to bridge the gap between complex theoretical models and practical engineering needs, enhancing the capability to predict and mitigate the adverse effects of thermal loads on spacecraft components and thus contributing to the reliability and performance of space missions.
This paper is structured as follows:
Section 2 is devoted to introducing the considered system and the addressed problem by also recalling basic concepts concerning the thermal behavior of orbiting systems and thermo-mechanical theory.
Section 3 describes the proposed simplified model at both kinematic and dynamic levels, while
Section 4 presents the results obtained via the fully coupled system and the comparison with a finite element formulation. The linearization of the system equations is also discussed, along with an approach to infer-thermal flutter. Finally,
Section 5 draws the work’s conclusions.
2. Problem Introduction: Thermo-Mechanical Modeling
A thermal shock associated with a sudden warming or cooling process is commonly accompanied, in a structure, by more or less pronounced dynamical effects. Depending on the geometry of the system and its thermal properties, such phenomena can be considered a negligible noise, a disturbance or a severe alteration of the expected performance. Among the many kinds of thermo-mechanical interaction, we shall consider in the following the well-known phenomenon of thermal jitter [
30], which usually occurs at the sun–shadow and shadow–sun transitions and affects light and slender structures that move along orbits with a period of eclipse.
A simple system was devised to highlight the basic mechanism of the interaction and used to draw conclusions from the proposed model. The spacecraft, located on a geostationary orbit (GEO), is in fact modeled as a thin plate of height
L, unit width and thickness
d, with
, as illustrated in
Figure 1. The
x-axis is normal to the plate and always Sun-pointing; the
y-axis is normal to the orbital plane. The scenario corresponds to an equinox when the eclipse duration is maximum, at about 70 min, and the temperature excursion has the maximum amplitude. The true anomaly
defines the instantaneous position of the spacecraft with respect to the subsolar point.
For a rough evaluation of the average temperature onboard, the spacecraft can be modeled as isothermal; a single thermal balance equation is, therefore, sufficient to plot the temperature history along the orbit.
Due to its simple shape, a unit area of the plate was considered. The subscripts
f and
b refer to the front and back of the plate with respect to the Sun,
m is the mass per unit length, and the other terms are the classical thermo-optical parameters. The lateral surface of the thickness was neglected.
Figure 2 shows the trend of the average temperature along the orbit, with rapid and pronounced cooling of the plate in the shadow zone due to its lightness. This significant change in temperature would result only in a length variation without altering the shape of the plate, and there would not be associated dynamical phenomena.
Nevertheless, it is not realistic to assume the same temperature for the surface facing the Sun and the opposite one, facing deep space. An improvement in the accuracy can be obtained by splitting the plate into two layers that are suitably connected. This permits to take into account the thermal gradient across the thickness d.
The mathematical model, in this case, consists of two thermal-balance equations coupled through the conductive heat exchange:
The equations are written for the unit area; the conductance is
, where
k is the thermal conductivity,
is the cross-section area of the honeycomb (in the case of a sandwich panel; for a homogeneous material, it is
). The specific heat,
c, is the same for the two layers.
The graphs in
Figure 3 and
Figure 4, corresponding to the improved modelization, highlight a non-negligible thermal gradient across the thickness. Two different combinations of the emissivities were considered. In Case A of
Figure 3, the low value of
reduces the outgoing heat flux so that the absorbed sun radiation increases the temperature of the front layer with respect to the back one. Notice that, during the cooling phase in the shadow zone, the gap reduces, but the temperature
always remains higher than
. When the emissivities are swapped, as in Case B of
Figure 4, most of the heat absorbed via the front layer is re-radiated and, therefore, lost. The thermal gradient is considerably reduced, and remarkably, the front layer, during the eclipse, cools faster than the back one.
This last aspect can be better highlighted by observing the temperature difference between the two faces, as illustrated in
Figure 5 and
Figure 6.
In Case A (
Figure 5), the temperature difference
, constant during the daytime, decreases very quickly at the beginning of the nighttime, in the cooling phase, and it tends asymptotically to zero at the end of the eclipse; it remains, however, always positive.
In Case B (
Figure 6), for the reasons stated above,
is much smaller than in the previous case, and it shows two overshoots in both the cooling and warming phases. Moreover, in the cooling phase, it becomes negative. It is interesting to notice that Case B represents a behavior very similar to that observed for the boom gradient temperatures of the Hubble Space Telescope, as presented in [
8], and it was considered the reference condition for this investigation.
In both cases, it is evident that, during the transit in the shadow zone and in the first portion of the orbit in sunlight, the panel undergoes a non-negligible thermal transient, and .
Thermostructural Coupling
To the variation of the thermal field is associated a change in geometry. For a monodimensional element, we have the basic relationship
In the first isothermal modeling, the thermal transients would produce only a contraction and elongation of the panel; the latter would not change its rectilinear shape. With the introduction of the thermal gradient associated with the refined model, the effect is different, and the panel assumes a curvature depending on
. Indeed, if we further improve the accuracy by increasing the number of layers along with the corresponding thermal-balance equations, we find a linear temperature distribution across the thickness. If we imagine the latter as a stack of frictionless layers, like a deck of cards, the temperature distribution would produce only a differential elongation, leaving them flat, as shown in
Figure 7a.
In a real, coherent material, the shear stresses
between adjacent layers cause a condition of compulsion; the hotter layer is inhibited in its elongation by the colder layer, whose elongation, corresponding to a lower temperature, is less. The consequent geometry shows a variation in shape, as illustrated in
Figure 7b. However, any intermediate layer is self-balanced; the
distribution has, in fact, a null resultant, as described in
Figure 7c. The plate, initially straight, acquires a curvature whose amount and sign are not constant since
. It is possible to find, for a given
, the expression of a thermal bending moment, equivalent to a mechanical moment, which produces the same curvature without the associated state of stress.
For this purpose, let us consider an element of the plate with the current temperature distribution. The latter can be considered as a combination of a constant distribution, the average temperature, and a butterfly distribution expressing the gradient, as depicted in
Figure 8.
From the above relationship,
where, in this case,
is the strain,
the stress and
E the Young modulus.
With reference to
Figure 8, the average temperature cannot modify the shape of the panel; if anything, only its length. The butterfly trend is instead responsible for the curvature. The expression for the linear distribution is
To obtain the expression of the thermal bending moment
, we consider the cross section of the plate,
Figure 9, where the normal stress is applied:
where
is just the moment of inertia of the cross section of the unit width. If
,
It should not be forgotten that, since
, also,
.
3. Simplified Model
In this section, the simplified model is introduced. The dynamics of the orbiting system is discretized and reduced to a mass–spring–dashpot paradigm. At first, the kinematic behavior of the system is discussed. Then, the fully coupled thermo-mechanical governing equations are presented.
3.1. Kinematics of the Panel
The exposure to solar radiation causes, therefore, the panel to move according to the forcing term . Such a kind of motion can be smartly modeled by means of modal analysis with boundary conditions depending on time. This approach, although elegant, has the drawback of being long and laborious. As an alternative, a simplified model is proposed below.
The continuous system with distributed parameters is discretized and reduced to a single degree of freedom with kinematics that resembles the first natural mode of vibration. The panel is replaced with two pieces of the same length, , connected with a cylindrical hinge. The origin O of the body-reference frame is located in the center of mass (c.m.) of the system. B and are located at the middle point of each piece.
Since there are no forces applied, the system keeps the momentum, and the c.m. cannot move; so, when the configuration changes due to the action of
, the point
A is constrained to move along the
x-axis and the points
B and
along the
y-axis, as illustrated in
Figure 10. This is the only way to satisfy the condition of the momentum.
The shape of the panel is governed according to ; when , and , , when , and , when , and again .
The kinematics of any other point of the segment can be described easily once a local abscissa
is introduced. Its origin is located in
A and assumed to be positive toward
B, as in
Figure 11. For any value of
, the generic point
Q of abscissa
will have coordinates in the body frame,
and
. These expressions can be derived starting from the vertical position of the element and applying first a rotation
around
A and then a translation along
x to relocate
B on the
y-axis as the constraint requires. When
,
and when
,
, in both cases, the expressions are coherent.
We observe that
and
are the parametric expressions of motion where
. Eliminating the dependence on time, we obtain the trajectory, followed by
Q. We can write
Moreover,
after some algebra, we get
When point
A, i.e., the hinge, moves back and forth along the
x-axis, any point of the element
, except for
B itself which moves along the
y-axis, describes an elliptical arc with semiaxes
and
. Due to the mirror symmetry with respect to the
x-axis, the same conclusions apply to the element
;
Figure 10.
The family of curves shown in
Figure 12 represents the trajectory of the point
Q (and its gemini on the element
), for
and
from zero to
m with a step of 2.5 m. Semimajor and semiminor axes swap with each other as
increases from zero to
. When
, the trajectory is the arc of a circle. These curves are mainly of academic interest since the expected values of
are much smaller.
At the end of this digression, the two-way coupling between thermal dynamics and structural dynamics is evident: the thermal gradient produces a thermal bending moment that curves the panel and modifies its projected area with respect to the Sun; this, in turn, affects the thermal gradient and so on. In addition, due to the low thermal inertia of the panel, the bending moment changes very rapidly, inducing inertial effects and triggering vibrations.
3.2. Dynamics of the Panel
The analysis of the phenomenon must, therefore, include the dynamics of the panel, together with the thermal balance equations.
We derived the dynamics equation following the Lagrangian approach. As a first step, we have to write the expression for the kinetic energy of the simplified model; it is, of course, a function of .
With reference to
Figure 11, the generic point
Q moves with a velocity whose components are
The kinetic energy associated with an infinitesimal element of length
is
where
m is the mass per unit length. The total energy of the element
is, then,
It yields (see also
Appendix A)
As a second step, to simulate the structural elasticity of the continuous panel, we introduce a concentrated spring of stiffness
, located in correspondence with the hinge. The potential energy, therefore, has the form
For completeness, we take into account the structural damping by introducing Rayleigh’s dissipation function,
where
is a constant damping coefficient.
The Lagrangian function is
, and the dynamics equation is provided by
Applying formalism, we get
For small values of
, the equation reproduces the damped harmonic motion.
4. The Fully Coupled System
The problem is then reformulated to highlight the interdependence between the different dynamics. The incoming solar flux is a function of the angle
(this term, of course, vanishes during the eclipse), and the forcing term
has been made explicit with respect to the variables
and
. The radiative exchange between the facing surfaces on the current, concave side of the panel is neglected.
The reduction in the degrees of freedom from a distributed-parameters model to a discrete one inevitably introduces a certain degree of arbitrariness in the choice of representative parameters, in particular concerning
, the stiffness of the concentrated spring in the hinge. Criteria can be identified to make the simplified model accurate in the description and simulation of specific physical aspects relevant to design engineering, such as the system attitude, the displacement excursion or the frequency of structural oscillations. Therefore, the selection of
is driven by the model user choices, according to the design needs. More details are also provided in
Section 4.3.
A criterion could be to compare the deflection in the midspan of the panels in the continuous and the discrete models, subjected to the same bending moment applied to the tips, as in
Figure 13. Other criteria are, nevertheless, eligible.
Given
, the deflection
in the midspan of the continuous panel is
In the discrete model, the deflection is
; imposing equality
and, therefore,
but, at the equilibrium,
, so that
By integrating Equation (
20) with a numerical solver, we get the plot of
Figure 14. We can see that, for most of the orbit, when the system is in thermal equilibrium, the moment
is constant, and so the effect it produces, the
angle is in fact constant with a value of about
. The thermo-mechanical properties of the system are listed in
Table 1.
Abruptly, at the entrance of the shadow zone, the thermal equilibrium is lost.
quickly reduces its magnitude and changes sign. The panel reverses its shape and starts to vibrate due to dynamic effects (please refer to
Figure 15). Even while still vibrating, with a thermal equilibrium—a new one—not yet achieved, a new transient begins at the exit from the shadows, as magnified in
Figure 16. The sign of
changes again and so the sign of
, with dynamic effects progressively dampened until the successive eclipse.
The curve of
Figure 17 refers to the case of no damping,
, with the thermal gradient influenced by the variation of the area exposed to the solar flux. The system shows thermal flutter with a periodic progressive increment of the amplitude of oscillation. If we remove the coupling between the thermal gradient and the angle
, assuming a constant projected area, the system behaves as shown in
Figure 18.
The oscillation triggered by the dynamical load due to the eclipse, in the absence of damping, maintains the residual amplitude attained at the end of the transient. In this case, the mechanism of thermal pumping cannot be established.
As a final note, we can observe that the jitter frequency comes from the natural frequency,
, of the linearized dynamics equation
where
is the moment of inertia of half of a panel computed with respect to its center (
B or
). The frequency is then
This suggests another criterion to select the arbitrary value of
, focusing on the system frequency behavior instead.
Indeed, to impose equality between
and the frequency of the first free–free natural mode of the continuous model, one can obtain
It can be noticed how the value of the stiffness is greater than the one obtained from the comparison between the static displacements in the midspan.
4.1. System Linearization
Until the magnitude of the parameters is chosen within a range of realistic values, the angular displacement remains small, and the usage of the linearized dynamics equation appears justified.
The linearization of the variable , nevertheless, cannot be extended to the thermal equations because the two-way coupling would be lost. The thermal gradient would result in a forcing term for the dynamics equation, depending only on the angle , and it would not be affected by the projected area, which is, instead, ruled by the angle .
When, therefore, the dependence of the incoming solar flux on the angle
is maintained, the model is reduced to the following equations:
The same detail of
Figure 15 and
Figure 16 is repeated in
Figure 19 and
Figure 20 (with higher magnification) by proposing a comparison between linear and non-linear dynamics. The curves are almost indistinguishable.
4.2. Thermal Flutter Prediction
It is worth noting that the possibility of thermal flutter could be inferred avoiding the direct integration of equations. We observe, in fact, that once the transient is extinct, the excursion of
is very small and, as a consequence, the expected variation of
is even smaller since this last depends on
. Also the thermal equations can therefore be linearized. We remind the reader that
with
a common linearization temperature for both
and
.
Assuming, for simplicity,
, and subtracting the second equation from the first one, we get
with the positions
where the units are
K/s and
1/s (note that we are considering a surface of a unit area that is transparent in the equation).
The linearized thermal equations become
and, with the additional grouping
the complete system can be written:
Due to the low thermal inertia, the response of the system, in terms of
, to the variation in the projected area, follows the dynamics of
just with a little delay. The first equation, then, can be written:
and the solution
The first term of the second member is the solution of the associated homogeneous equations; it takes into account the initial condition but vanishes over time and is not interesting. The second term is the long-period solution.
After some time, the only surviving terms are
The forcing term of the dynamics has the same frequency as the dynamics itself; the system is resonant. Equation (
34) has a solution in closed form, and the amplitude of the oscillation can, however, be controlled with a suitable rate of damping.
4.3. Comparison with Finite-Element Approach
The results obtained via the formulation presented in Equation (
28) are here compared and validated with a model of the orbiting structure realized using a finite-element (FEM) approach. The structure, illustrated in
Figure 21, is a beam realized in aluminum material in an MSC Nastran suite, having dimensions and properties equal to those listed in
Table 1. The system is discretized via Bar2 elements for a total of 101 nodes. A thermal bending moment,
, is obtained from the full system dynamics simulations and imported into the FEM environment. A modal transient response simulation is performed by applying the time-variant moment to the structure in order to retrieve the resulting nodal displacements and rotations.
To reproduce the in-orbit free–free dynamics, subjected to the thermal bending moment
, and compute at the same time structural displacements, the inertia relief (IR) [
31] technique was adopted in the FE model. Indeed, conventional static analysis cannot be carried out for unconstrained systems, as the stiffness matrix would be singular. On the other hand, the IR method utilizes a different approach to obtain the load equilibrium in the model, with no applied constraints. IR includes inertial forces experienced by the body in the solution through inertia-relief loading that balances the external loading in a three-step procedure: (i) rigid body-mode evaluation, (ii) the assessment of inertia-relief loads, and (iii) the solution of a supported, self-equilibrated static load case. In particular, the IR is applied at the beam’s center of mass in node
, as represented in
Figure 21.
To compare the results of FEM with the simplified model, a
angle can be computed according to the following relation:
where
is the total deflection of the beam, obtained via the algebraic sum of the tip and center nodes’ displacements. Such an angle can be directly compared with the variable
obtained from the simplified model.
Some considerations can be drawn regarding the selection of the elastic stiffness
in the simplified model.
Figure 22 illustrates a comparison between
and
when
is determined according to the deflection criterion outlined in Equation (
24). The two approaches show good agreement, as indicated by the matching deflections.
On a different note, it can be observed—as expected—that the frequency of the jitter phenomenon is not perfectly reproduced by selecting the stiffness
to match the maximum deflection of the slender structure. Indeed, the first natural frequency of the beam can be computed according to Equation (
27)—or via the FEM model—as
Hz. As the thermal shock acts on the structure by exciting its first bending mode, one should expect the jitter to present mainly this frequency content. As can be noticed in
Figure 23, the frequency of the jitter shows a mismatch of about 12% with respect to
. On the other hand, if the designer selects
, according to the criterion in Equation (
27), the frequency content aligns perfectly with the expected mode (see
Figure 24), bearing the cost of a discrepancy between
and
of a 0.93-deg RMSE (root mean square error).
In both the considered cases, however, is proportional to the flexural rigidity, , of the beam, with the difference of a factor . Consequently, if selecting a in the simplified model, it is reasonable to assume that the flexural stiffness of the structure is increased—with the length being the same. This can be reflected also in the thermo-mechanical coupling term (the forcing term on the one-d.o.f. model) by increasing the term D by the factor . After this operation, one can have a model that provides, at the same time, an angle response and a frequency content aligned with the FEM model results (both the maximum deflection and frequency criteria are met). It is, however, left to the user which aspect to prioritize in the definition of the model parameters.
5. Conclusions
This paper has developed a simplified thermo-mechanical model of a two-layer orbiting slender structure for thermal jitter and flutter evaluations. A one-degree-of-freedom model, defined according to a mass–spring–damper paradigm, is presented to replicate the dominant impact of the thermal shock and jitter phenomena on the attitude of the system, according to maximum structural deflection and the first natural mode of vibration. Based on the proposed model, the temperature difference across the structure was calculated, and thermal vibrations induced via the change from sunlight to umbra conditions and vice versa were evaluated for a flexible spacecraft in a GEO environment. Some criteria and guidelines on how to select the simplified model parameters according to physical considerations are provided for engineering applications while also comparing the model outcome with an equivalent approach via the finite element method. The results reveal that the model could be useful for forecasting thermally induced vibrations and flutter conditions for aerospace applications and for understanding the energy pumping mechanism that can affect space missions, leading to attitude instability.
The presented approach not only addresses a gap in the available thermal design instruments but also offers an accessible solution for tackling thermal jitter problems in spacecraft design. As such, it represents an additional approach to be used in the field of thermal–structural dynamics, with the potential for further development and broader application. In future research, the approach could be extended to diverse structural configurations, including more layers in the model to accurately reproduce composite materials, broadening its relevance to a wider range of spacecraft designs. Additionally, incorporating adaptive modeling techniques or scheduling could allow for dynamic adjustments of the model parameters, improving accuracy in varying thermal environments. Further validation through comparison with more complex simulations—such as multi-physics thermo-structural simulations to verify thermal flutter occurrences—and, eventually, experimental validation would also strengthen confidence in the method’s robustness and reliability. Ultimately, this approach has the potential to be integrated into comprehensive spacecraft design tools, offering engineers a quick and efficient means of predicting and mitigating thermal flutter and jitter effects across various mission scenarios.