Next Article in Journal
Validation of Experimental Data for the Application of the Magnesium Alloy “Elektron 43”
Previous Article in Journal
A Physical and Spectroscopic Survey of the Lunar South Pole with the Galileo Telescope of the Asiago Astrophysical Observatory
Previous Article in Special Issue
Enabling High-Power Conditioning and High-Voltage Bus Integration Using Series-Connected DC Transformers in Spacecrafts
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Thermo-Mechanical Jitter in Slender Space Structures: A Simplified Modeling Approach

by
Maurizio Parisse
and
Federica Angeletti
*
School of Aerospace Engineering, University of Rome “La Sapienza”, 00138 Rome, Italy
*
Author to whom correspondence should be addressed.
Aerospace 2024, 11(9), 694; https://doi.org/10.3390/aerospace11090694
Submission received: 8 July 2024 / Revised: 22 August 2024 / Accepted: 23 August 2024 / Published: 25 August 2024
(This article belongs to the Special Issue Advanced Spacecraft/Satellite Technologies)

Abstract

:
Thermally induced vibrations usually affect spacecraft equipped with light and slender appendages such as booms, antennas or solar panels. This phenomenon occurs when a thermal shock, resulting from the sudden cooling and warming phases at the entrance and exit from eclipses, triggers mechanical vibrations. The study proposed hereafter concerns the modeling and prediction of jitter of thermal origin in a long and thin plate with a sun-pointing attitude in geostationary orbit. The system’s temperature and dynamics are described by a set of equations expressing the two-way coupling between the thermal bending moment and the shape of the panel. The structure is discretized and reduced to a one-degree-of-freedom simplified model able to identify a mechanism of thermal pumping that could lead to instability. Finally, the results of the analysis are compared with those obtained with a more accurate FEM modelization.

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 L > > 1 > > d , 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.
m c d T d t = α C σ ( ε f + ε b ) T 4 .
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:
1 2 m c d T f d t = α C K ( T f T b ) σ ε f T f 4 1 2 m c d T b d t = K ( T f T b ) σ ε b T b 4 .
The equations are written for the unit area; the conductance is K = k A c / d , where k is the thermal conductivity, A c is the cross-section area of the honeycomb (in the case of a sandwich panel; for a homogeneous material, it is A c = A = 1 ). 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 ε f 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 T f always remains higher than T b . 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 Δ T = T f T b , 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, Δ T 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 Δ T = Δ T ( t ) .

Thermostructural Coupling

To the variation of the thermal field is associated a change in geometry. For a monodimensional element, we have the basic relationship
Δ L = L L 0 = a L 0 ( T T 0 ) .
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 Δ T . 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 Δ T = Δ T ( t ) . It is possible to find, for a given Δ T , 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,
Δ L L 0 = L L 0 L 0 = ε = a Δ T σ = ε E = a E Δ T
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
T ( x ) = T f T b d x .
To obtain the expression of the thermal bending moment M T , we consider the cross section of the plate, Figure 9, where the normal stress is applied:
σ = a E T f T b d x
M T = d 2 d 2 σ x d x = a E T f T b d d 2 d 2 x 2 d x = a E T f T b d d 3 12
where d 3 / 12 = I is just the moment of inertia of the cross section of the unit width. If ( T f T b ) / d = G ,
M T = a E G I .
It should not be forgotten that, since Δ T = T f T b = Δ T ( t ) , also, M T = M T ( t ) .

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 M T = M T ( t ) . 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, L / 2 = 2 l , 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 B 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 M T ( t ) , the point A is constrained to move along the x-axis and the points B and B 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 ϑ = 0 , O A and B ( 0 , l ) , B ( 0 , l ) , when ϑ = 90 , A ( l , 0 ) and O B B , when ϑ = 90 , A ( l , 0 ) and again O B B .
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 ϑ > 0 , the generic point Q of abscissa ζ will have coordinates in the body frame, x = ( l ζ ) sin ϑ and y = ζ cos ϑ . 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 ζ = 0 , Q A and when ζ = l , Q B , in both cases, the expressions are coherent.
We observe that x = ( l ζ ) sin ϑ and y = ζ cos ϑ are the parametric expressions of motion where ϑ = ϑ ( t ) . Eliminating the dependence on time, we obtain the trajectory, followed by Q. We can write
y = ζ cos ϑ = ζ 1 sin 2 ϑ but sin 2 ϑ = x 2 ( l ζ ) 2
Moreover,
y 2 = ζ 2 ( 1 sin 2 ϑ ) = ζ 2 1 x 2 ( l ζ ) 2
after some algebra, we get
x 2 ( l ζ ) 2 + y 2 ζ 2 = 1 .
When point A, i.e., the hinge, moves back and forth along the x-axis, any point of the element A B , except for B itself which moves along the y-axis, describes an elliptical arc with semiaxes l ζ and ζ . Due to the mirror symmetry with respect to the x-axis, the same conclusions apply to the element A B ; Figure 10.
The family of curves shown in Figure 12 represents the trajectory of the point Q (and its gemini on the element A B ), for 90 ϑ 90 and ζ from zero to 2 l = 25 m with a step of 2.5 m. Semimajor and semiminor axes swap with each other as ζ increases from zero to 2 l . When ζ = l / 2 , 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 M T 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 ϑ ( t ) .
With reference to Figure 11, the generic point Q moves with a velocity whose components are
x ˙ = ( l ζ ) ϑ ˙ cos ϑ and y ˙ = ζ ϑ ˙ sin ϑ .
The kinetic energy associated with an infinitesimal element of length d ζ is
d E k = 1 2 m ( x ˙ 2 + y ˙ 2 ) d ζ
where m is the mass per unit length. The total energy of the element A B is, then,
E k = 0 2 l 1 2 m ( x ˙ 2 + y ˙ 2 ) d ζ = 1 2 m 0 2 l [ ( l ζ ) 2 ϑ ˙ 2 cos 2 ϑ + ζ 2 ϑ ˙ 2 sin 2 ϑ ] d ζ
It yields (see also Appendix A)
E k = m l 3 3 ϑ ˙ 2 ( 1 + 3 sin 2 ϑ ) .
As a second step, to simulate the structural elasticity of the continuous panel, we introduce a concentrated spring of stiffness K e l , located in correspondence with the hinge. The potential energy, therefore, has the form
U = 1 2 K e l ϑ 2 .
For completeness, we take into account the structural damping by introducing Rayleigh’s dissipation function,
F = 1 2 c d ϑ ˙ 2
where c d is a constant damping coefficient.
The Lagrangian function is L = E k U , and the dynamics equation is provided by
d d t L ϑ ˙ L ϑ + F ϑ ˙ = M T
Applying formalism, we get
ϑ ¨ = 1 ( 1 + 3 sin 2 ϑ ) 3 M T 2 m l 3 3 ϑ ˙ 2 sin ϑ cos ϑ 3 c d 2 m l 3 ϑ ˙ 3 K e l 2 m l 3 ϑ .
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 M T has been made explicit with respect to the variables T f and T b . The radiative exchange between the facing surfaces on the current, concave side of the panel is neglected.
T ˙ f = 2 m c α C cos ϑ K ( T f T b ) σ ε f T f 4 T ˙ b = 2 m c K ( T f T b ) σ ε b T b 4 ϑ ¨ = 1 ( 1 + 3 sin 2 ϑ ) 3 a E I 2 m l 3 d ( T f T b ) 3 ϑ ˙ 2 sin ϑ cos ϑ 3 c d 2 m l 3 ϑ ˙ 3 K e l 2 m l 3 ϑ .
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 K e l , 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 K e l 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 M T , the deflection η in the midspan of the continuous panel is
η = M T L 2 8 E I .
In the discrete model, the deflection is 2 l sin ϑ ; imposing equality
2 l sin ϑ = M T L 2 8 E I
and, therefore,
ϑ = arcsin M T L 4 E I
but, at the equilibrium, M T = K e l ϑ , so that
K e l = M T arcsin M T L 4 E I if ϑ < < 1 K e l = 4 E I L .
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 M T is constant, and so the effect it produces, the ϑ angle is in fact constant with a value of about 5 . 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. M T 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 M T 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, c d = 0 , 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, ω n , of the linearized dynamics equation
I B ϑ ¨ + c d ϑ ˙ + K e l ϑ = M T
where I B = 2 m l 3 / 3 is the moment of inertia of half of a panel computed with respect to its center (B or B ). The frequency is then
ω n = K e l I B .
This suggests another criterion to select the arbitrary value of K e l , focusing on the system frequency behavior instead.
Indeed, to impose equality between ω n and the frequency of the first free–free natural mode of the continuous model, one can obtain
K e l I B = ( 1.506 π ) 2 E I m L 4 K e l = ( 1.506 π ) 4 E I m L 4 I B 5.21 E I L
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:
T ˙ f = 2 m c α C cos ϑ K ( T f T b ) σ ε f T f 4 T ˙ b = 2 m c K ( T f T b ) σ ε b T b 4 ϑ ¨ = 1 I B a E I d ( T f T b ) c d ϑ ˙ K e l ϑ .
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 Δ T is even smaller since this last depends on cos ϑ . Also the thermal equations can therefore be linearized. We remind the reader that
T 4 = ( T * + T T * ) 4 = T * 4 + 4 T * 3 ( T T * ) + 4 T * 3 T 3 T * 4
with T * a common linearization temperature for both T f and T b .
Assuming, for simplicity, ε f = ε b = ε , and subtracting the second equation from the first one, we get
T ˙ f T ˙ b = 2 α C m c cos ϑ 4 K + 8 σ ε T * 3 m c ( T f T b )
with the positions
Δ T ˙ = T ˙ f T ˙ b p = 2 α C m c q = 4 K + 8 σ ε T * 3 m c
where the units are [ p ] = K/s and [ q ] = 1/s (note that we are considering a surface of a unit area that is transparent in the equation).
The linearized thermal equations become
Δ T ˙ + q Δ T = p cos ϑ
and, with the additional grouping
r = c d I B s = a E I I B d
the complete system can be written:
Δ T ˙ + q Δ T = p cos ϑ ϑ ¨ + r ϑ ˙ + ω n 2 ϑ = s Δ T
Due to the low thermal inertia, the response of the system, in terms of Δ T , to the variation in the projected area, follows the dynamics of ϑ just with a little delay. The first equation, then, can be written:
Δ T ˙ + q Δ T = p cos ( ω n t + φ )
and the solution
Δ T ( t ) = e ( t t 0 ) q Δ T ( t 0 ) + t 0 t e ( t τ ) q p cos ( ω n τ + φ ) d τ
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.
Δ T ( t ) = p q 2 q 2 + ω n 2 { 1 q cos ( ω n t + φ ) e q t cos φ + ω n q 2 sin ( ω n t + φ ) e q t sin φ
After some time, the only surviving terms are
Δ T ( t ) = p q 2 q 2 + ω n 2 1 q cos ( ω n t + φ ) + ω n q 2 sin ( ω n t + φ )
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, M T , 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 M T , 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 O S , as represented in Figure 21.
To compare the results of FEM with the simplified model, a ϑ f angle can be computed according to the following relation:
ϑ f = arcsin ( 2 L η f )
where η f 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 K e l in the simplified model. Figure 22 illustrates a comparison between ϑ and ϑ f when K e l 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 K e l 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 f n = 0.144 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 f n . On the other hand, if the designer selects K e l , 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 ϑ ( t ) and ϑ f ( t ) of a 0.93-deg RMSE (root mean square error).
In both the considered cases, however, K e l is proportional to the flexural rigidity, D = E I , of the beam, with the difference of a factor χ = 5.21 4 = 1.3 . Consequently, if selecting a K e l = 5.21 D / L 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 M T (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.

Author Contributions

Conceptualization, M.P. and F.A.; Methodology, M.P. and F.A.; Software, M.P. and F.A.; Validation, M.P. and F.A.; Formal analysis, M.P. and F.A.; Investigation, M.P. and F.A.; Writing—original draft, M.P. and F.A.; Writing—review & editing, M.P. and F.A.; Supervision, M.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Data is contained within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Appendix A

The expression of the kinetic energy can be derived directly, avoiding the differential approach. The motion of the element A B is, in fact, a combination of a rotation around B, together with a rigid translation along the y-axis. To respect the constraint, the angular velocity, ϑ ˙ , and the velocity of translation, y ˙ B , are not independent.
The kinetic energy is, therefore, the sum of the translational and rotational kinetic energy: E k = E T + E R .
E T = 1 2 2 l m y ˙ B 2 = m l ( l ϑ ˙ sin ϑ ) 2 = m l 3 ϑ ˙ 2 sin 2 ϑ E R = 1 2 I B ϑ ˙ 2 = m l 3 3 ϑ ˙ 2
where I B = 2 m l 3 / 3 is the moment of inertia of the element A B with respect to the axis passing through B and normal to the plane of motion. We have
E k = E T + E R = m l 3 3 ϑ ˙ 2 ( 1 + 3 sin 2 ϑ ) .

References

  1. Quadrelli, M.B.; Mettler, E.; Soloway, D.; Kelkar, A. Controls structure interaction on the Jupiter Icy Moons Orbiter. Acta Astronaut. 2009, 65, 766–786. [Google Scholar] [CrossRef]
  2. Angeletti, F.; Tortorici, D.; Laurenzi, S.; Gasbarri, P. Vibration Control of Innovative Lightweight Thermoplastic Composite Material via Smart Actuators for Aerospace Applications. Appl. Sci. 2023, 13, 9715. [Google Scholar] [CrossRef]
  3. Angeletti, F.; Iannelli, P.; Gasbarri, P.; Gonzalez, J.A.P.; Ellero, N.; Wattrelot, T.; Ankersen, F.; Sabatini, M.; Celani, F.; Palmerini, G. Robust Collocated Control of Large Flexible Space Structures. IFAC—PapersOnLine 2022, 55, 85–90. [Google Scholar] [CrossRef]
  4. Oda, M.; Hagiwara, Y.; Suzuki, S.; Nakamura, T.; Inaba, N.; Sawada, H.; Yoshii, M.; Goto, N. Measurement of Satellite Solar Array Panel Vibrations Caused by Thermal Snap and Gas Jet Thruster Firing. In Recent Advances in Vibrations Analysis; Baddour, N., Ed.; IntechOpen: Rijeka, Croatia, 2011; Chapter 7. [Google Scholar] [CrossRef]
  5. Etkin, B.; Hughes, P.C. Explanation of the Anomalous Spin Behaviour of Satellites with Long Flexible Antenna. J. Spacecr. Rocket. 1967, 4, 1139–1145. [Google Scholar] [CrossRef]
  6. Kershner, R.B.; Fischell, R.E. Gravity-Gradient Stabilization of Earth Satellites. IFAC Proc. Vol. 1965, 2, 249–266. [Google Scholar] [CrossRef]
  7. Denney, C.; Welch, R.; Zimbelman, D. Sunrise/sunset thermal shock disturbance analysis and simulation forthe TOPEX satellite. In Proceedings of the 28th Aerospace Sciences Meeting, Reno, NV, USA, 8–11 January 1990. [Google Scholar] [CrossRef]
  8. Foster, C.; Tinker, M.; Nurre, G.; Till, W. The Solar Array-Induced Disturbance of the Hubble Telescope Pointing System; NASA Technical Paper 3556; MSFC, NASA: Huntsville, AL, USA, 1995. [Google Scholar]
  9. Boley, B.A. Thermally Induced Vibrations of Beams. J. Aeronaut. Sci. 1956, 23, 179–181. [Google Scholar]
  10. Boley, B.A.; Barber, A.D. Dynamic Response of Beams and Plates to Rapid Heating. J. Appl. Mech. 1957, 24, 413–416. [Google Scholar] [CrossRef]
  11. Boley, B.A. Approximate Analyses of Thermally Induced Vibrations of Beams and Plates. J. Appl. Mech. 1972, 39, 212–216. [Google Scholar] [CrossRef]
  12. Thornton, E.A.; Kim, Y.A. Thermally induced bending vibrations of a flexible rolled-up solar array. J. Spacecr. Rocket. 1993, 30, 438–448. [Google Scholar] [CrossRef]
  13. Frisch, H.P. Thermally induced vibrations of long thin-walled cylinders of open section. J. Spacecr. Rocket. 1970, 7, 897–905. [Google Scholar] [CrossRef]
  14. Tsuchiya, K. Thermally induced vibrations of a flexible appendage attached to a spacecraft. AIAA J. 1977, 15, 505–510. [Google Scholar] [CrossRef]
  15. Li, J.; Yan, S.; Cai, R. Thermal analysis of composite solar array subjected to space heat flux. Aerosp. Sci. Technol. 2013, 27, 84–94. [Google Scholar] [CrossRef]
  16. Johnston, J.D.; Thornton, E.A. Thermally Induced Dynamics of Satellite Solar Panels. J. Spacecr. Rocket. 2000, 37, 604–613. [Google Scholar] [CrossRef]
  17. Su, X.; Zhang, J.; Wang, J.; Bi, Y.; Qie, D.; Xiang, Z.; Xue, M. Experimental investigation of the thermally induced vibration of a space boom section. Sci. China Phys. Mech. Astron. 2015, 58, 1–9. [Google Scholar] [CrossRef]
  18. Mason, J.B. Analysis of Thermally Induced Structural Vibrations by Finite Element Techniques; NASA Technical Report; GSFC, NASA: Greenbelt, MD, USA, 1968. [Google Scholar]
  19. Chen, X.; Mohan, R.V.; Tamma, K.K. Instantaneous response of elastic thin-walled structures to rapid heating. Int. J. Numer. Methods Eng. 1994, 37, 2389–2408. [Google Scholar] [CrossRef]
  20. Appel, S.; Wijker, J. Structural Modelling for Thermoelastic Analysis. In Simulation of Thermoelastic Behaviour of Spacecraft Structures: Fundamentals and Recommendations; Springer International Publishing: Cham, Switzerland, 2022; pp. 111–163. [Google Scholar] [CrossRef]
  21. Xue, M.; Ding, Y. Two kinds of tube elements for transient thermal–structural analysis of large space structures. Int. J. Numer. Methods Eng. 2004, 59, 1335–1353. [Google Scholar] [CrossRef]
  22. Wang, J.; Jin, D.; Fan, C.; Bi, Y.; Su, X.; Liu, G.; Xiang, Z. Predicting the on-orbit thermally induced vibration through the integrated numerical and experimental approach. Acta Astronaut. 2022, 192, 341–350. [Google Scholar] [CrossRef]
  23. Shen, Z.; Hu, G. Thermally induced vibrations of solar panel and their coupling with satellite. Int. J. Appl. Mech. 2013, 05, 1350031. [Google Scholar] [CrossRef]
  24. Shen, Z.; Li, H.; Liu, X.; Hu, G. Thermal shock induced dynamics of a spacecraft with a flexible deploying boom. Acta Astronaut. 2017, 141, 123–131. [Google Scholar] [CrossRef]
  25. Mu, R.; Tan, S.; Wu, Z.; Qi, Z. Coupling dynamics of super large space structures in the presence of environmental disturbances. Acta Astronaut. 2018, 148, 385–395. [Google Scholar] [CrossRef]
  26. Cao, Y.; Cao, D.; He, G.; Liu, L. Thermal alternation induced vibration analysis of spacecraft with lateral solar arrays in orbit. Appl. Math. Model. 2020, 86, 166–184. [Google Scholar] [CrossRef]
  27. Liu, L.; Sun, S.; Cao, D.; Liu, X. Thermal-structural analysis for flexible spacecraft with single or double solar panels: A comparison study. Acta Astronaut. 2019, 154, 33–43. [Google Scholar] [CrossRef]
  28. Ganilova, O.A.; Cartmell, M.P.; Kiley, A. Application of a dynamic thermoelastic coupled model for an aerospace aluminium composite panel. Compos. Struct. 2022, 288, 115423. [Google Scholar] [CrossRef]
  29. Ganilova, O.A.; Cartmell, M.P.; Kiley, A. Experimental investigation of the thermoelastic performance of an aerospace aluminium honeycomb composite panel. Compos. Struct. 2021, 257, 113159. [Google Scholar] [CrossRef]
  30. Thornton, E.A. Thermal Structures for Aerospace Applications; AIAA Education Series; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 1996. [Google Scholar]
  31. Barnett, A.R.; Widrick, T.W.; Ludwiczak, D.R. Closed-form Static Analysis with Inertia Relief and Displacement-Dependent Loads Using a MSC/NASTRAN DMAP Alter. In Proceedings of the 1995 World Users’ Conference, Los Angeles, CA, USA, 8–12 May 1995. [Google Scholar]
Figure 1. Case study: slender plate with unit width, height L and thickness d orbiting the Earth, maintaining a Sun-pointing attitude.
Figure 1. Case study: slender plate with unit width, height L and thickness d orbiting the Earth, maintaining a Sun-pointing attitude.
Aerospace 11 00694 g001
Figure 2. Time history of average temperature.
Figure 2. Time history of average temperature.
Aerospace 11 00694 g002
Figure 3. Case A: temperature–time history for ε f < ε b .
Figure 3. Case A: temperature–time history for ε f < ε b .
Aerospace 11 00694 g003
Figure 4. Case B: temperature–time history for ε f > ε b .
Figure 4. Case B: temperature–time history for ε f > ε b .
Aerospace 11 00694 g004
Figure 5. Magnification: temperature difference ( ε f < ε b ).
Figure 5. Magnification: temperature difference ( ε f < ε b ).
Aerospace 11 00694 g005
Figure 6. Magnification: temperature difference ( ε f > ε b ).
Figure 6. Magnification: temperature difference ( ε f > ε b ).
Aerospace 11 00694 g006
Figure 7. Thermal bending and stresses: (a) ideal frictionless layers, (b) shear stress effect on layers, (c) self-balancing.
Figure 7. Thermal bending and stresses: (a) ideal frictionless layers, (b) shear stress effect on layers, (c) self-balancing.
Aerospace 11 00694 g007
Figure 8. Temperature distribution across the structure.
Figure 8. Temperature distribution across the structure.
Aerospace 11 00694 g008
Figure 9. Normal stress in the structure cross section.
Figure 9. Normal stress in the structure cross section.
Aerospace 11 00694 g009
Figure 10. Simplified model: two parts of length 2 l connected with a cylindrical hinge.
Figure 10. Simplified model: two parts of length 2 l connected with a cylindrical hinge.
Aerospace 11 00694 g010
Figure 11. Motion components: projections on the coordinate axes.
Figure 11. Motion components: projections on the coordinate axes.
Aerospace 11 00694 g011
Figure 12. Trajectories’ envelope.
Figure 12. Trajectories’ envelope.
Aerospace 11 00694 g012
Figure 13. Panel deflection: continuous model (left) and discrete model (right).
Figure 13. Panel deflection: continuous model (left) and discrete model (right).
Aerospace 11 00694 g013
Figure 14. Thermal shock of the panel: angle ϑ .
Figure 14. Thermal shock of the panel: angle ϑ .
Aerospace 11 00694 g014
Figure 15. Magnification: entrance in the shadow zone.
Figure 15. Magnification: entrance in the shadow zone.
Aerospace 11 00694 g015
Figure 16. Magnification: exit from the shadow zone.
Figure 16. Magnification: exit from the shadow zone.
Aerospace 11 00694 g016
Figure 17. Thermal flutter: projected area depending on ϑ .
Figure 17. Thermal flutter: projected area depending on ϑ .
Aerospace 11 00694 g017
Figure 18. Thermal flutter: constant projected area.
Figure 18. Thermal flutter: constant projected area.
Aerospace 11 00694 g018
Figure 19. Entrance in the shadow zone: linearized vs. nonlinear.
Figure 19. Entrance in the shadow zone: linearized vs. nonlinear.
Aerospace 11 00694 g019
Figure 20. Exit from the shadow zone: linearized vs. nonlinear.
Figure 20. Exit from the shadow zone: linearized vs. nonlinear.
Aerospace 11 00694 g020
Figure 21. Finite-element free–free beam model subjected to thermal moment M T .
Figure 21. Finite-element free–free beam model subjected to thermal moment M T .
Aerospace 11 00694 g021
Figure 22. Comparison between FEM and simplified model approach ( K e l = 4 E I / L ).
Figure 22. Comparison between FEM and simplified model approach ( K e l = 4 E I / L ).
Aerospace 11 00694 g022
Figure 23. FFT of ϑ angle time variation ( K e l = 4 E I / L ).
Figure 23. FFT of ϑ angle time variation ( K e l = 4 E I / L ).
Aerospace 11 00694 g023
Figure 24. FFT of ϑ angle–time variation ( K e l = 5.21 E I / L ).
Figure 24. FFT of ϑ angle–time variation ( K e l = 5.21 E I / L ).
Aerospace 11 00694 g024
Table 1. Thermo-mechanical properties of the orbiting structure.
Table 1. Thermo-mechanical properties of the orbiting structure.
VariableSymbolValueUnits of Measure
Emissivity—front panel ε f 0.8-
Emissivity—back panel ε b 0.3-
Absorptivity α 0.8-
Specific heatc900J/kg K
Linear densitym5.51kg/m
Thermal expansion coefficienta 2.5 × 10 5 1/K
Total beam lengthL50m
Quarter of beam lengthl12.5m
Young modulusE70GPa
Thicknessd 0.04 m
Widthw1m
ConductanceK25W/K
Elastic stiffness K e l 4438Nm/rad
Bending inertiaI 8 × 10 7 m4
Damping c d 3Nms/rad
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Parisse, M.; Angeletti, F. Thermo-Mechanical Jitter in Slender Space Structures: A Simplified Modeling Approach. Aerospace 2024, 11, 694. https://doi.org/10.3390/aerospace11090694

AMA Style

Parisse M, Angeletti F. Thermo-Mechanical Jitter in Slender Space Structures: A Simplified Modeling Approach. Aerospace. 2024; 11(9):694. https://doi.org/10.3390/aerospace11090694

Chicago/Turabian Style

Parisse, Maurizio, and Federica Angeletti. 2024. "Thermo-Mechanical Jitter in Slender Space Structures: A Simplified Modeling Approach" Aerospace 11, no. 9: 694. https://doi.org/10.3390/aerospace11090694

APA Style

Parisse, M., & Angeletti, F. (2024). Thermo-Mechanical Jitter in Slender Space Structures: A Simplified Modeling Approach. Aerospace, 11(9), 694. https://doi.org/10.3390/aerospace11090694

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop