Xiaode YUAN,Zhihai XIANG

Department of Engineering Mechanics,School of Aerospace Engineering,Tsinghua University,Beijing 100084,China

KEYWORDS Fourier finite element;Stability criteria;Thermal flutter;Thermally Induced Vibration(TIV);Thin-walled structures

Abstract The flexible attachments of spacecraft may undergo Thermally Induced Vibration(TIV)on orbit due to the suddenly changed solar heating.The unstable TIV,called thermal- flutter,can cause serious damage to the spacecraft.In this paper,the coupled bending-torsion thermal vibration equations for an open thin-walled circular cantilever beam are established.By analyzing the stability of these equations based on the first Lyapunov method,the thermal- flutter criterion can be obtained.The criterion is very different form that of closed thin-walled beams because the torsion has great impact on the stability of the TIV for open thin-walled beams.Several numerical simulations are conducted to demonstrate that the theoretical predictions agree very well with the finite element results,which mean that the criterion are reliable.

1.Introduction

The flexible attachments of spacecraft generally have the characteristics of large size,light weight,low stiffness and small heat capacity.Therefore,these structures are prone to experiencing the Thermally Induced Vibration(TIV)due to the suddenly applied solar heat flux when the spacecraft enter or leave eclipse.1–4These vibrations could reduce the pointing accuracy of spacecraft and even introduce damage into the structure,especially when the vibration is unstable,i.e.,the thermal flutter.

TIV was firstly predicted theoretically by Boley as early as 1956.5Boley and Barber6showed that when a very thin beam or plate is subjected to rapid surface heating,the vibration can be induced by a kind of time-dependent thermal moment due to the rapid temperature gradient in the structure.Later on,the Boley parameter B= τTω1was defined to characterize the severity of TIV for cantilever beams,7where τTis the thermal characteristic time and ω1is the minimum angular frequency of the beam.The ratio of the maximum dynamic deflection over the quasi-static deflection of a cantilever beam can be expressed aswhich means that the smaller B is,the more severe the TIV is.Although the Boley parameter B is a nice index for pure bending TIV of a cantilever beam,practical structures may have more complex TIV modes.For example,the structure composed of open thin-walled beams is apt to undergo torsional vibration due to its ultra-low torsional stiffness.8

Compared to stable TIV,the thermal flutter is more harmful to space structures.This phenomenon was first observed on orbit in 19689and then it was realized in a laboratory environment.10After that,more and more coupled thermal-structure analyses were conducted to investigate the condition of thermal flutter.11,12Yu first established the stability criterion on the TIV of a closed thin-walled cantilever beam subject to solar heating13and then that criterion was updated by Graham.14In Graham’s criterion,the thermal flutter will only happen when the beam axis points away from the sun,where the beam axis is defined as the vector pointing from the fixed end of the beam to the free end of the beam.An important conclusion of this criterion is that the normal-incident heat flux will not induce thermal flutter,which is contradictory to both experiment results15and numerical simulations.16Realizing that the stability analysis should be established on the deformed steady state instead of the original configuration of the beam,Zhang and Xiang proposed a new criterion,which conforms with the experimental and numerical results.17

All existing criteria of thermal flutter are only applicable to closed thin-walled beams.In contrast,the criterion for open thin-walled beams must consider the bending and torsion coupling deformations.Consequently,the circumferential incident angle of the heat flux should have great impact on the stability of the TIV.With a full consideration of these two points,this paper establishes a thermal- flutter criterion suitable for open thin-walled circular cantilever beams based on the first Lyapunov method.18

2.Coupled thermal-structural dynamic analysis

2.1.Analysis model and basic assumptions

As Fig.1 shows,two sets of coordinate systems are defined to describe the deformation of the cantilever beam.OXYZ is a fixed spatial coordinate system,in which X axis is the initial centroid axis and Y axis points to the initial opening direction.Oxyz is a local coordinate system attached at a point on the beam,in which x axis is the deformed centroid axis and y axis always points to the opening direction of the rotated beam.

The dimensions of the interested beam are defined as follows:l is the beam length;R and h are the midline radius and thickness of the beam cross-section,respectively.For a thin-walled slender beam,h/R?1 and R/l?1,so that Euler-Bernoulli beam theory is applicable.

The solar heat flux vector S0is uniformly distributed along the beam length.θ0is the angle between S0and vector n,which is the normal of the beam and opposite to the projection of S0in plane OYZ.α is the angle between S0and the Y axis in the YOZ plane.

The following assumptions are adopted in the analysis:

Fig.1 An open thin-walled circular cantilever beam subject to solar heat flux.

(1)Emission and radiation of heat to space is considered but convection and radiation between the different surfaces of the beam are neglected.

(2)Heat transfer along beam length is neglected.

(3)At positions of X=0,X=l and the longitudinal opening sides of the beam are adiabatic.

(4)The amplitude of the perturbation temperature is much smaller than the average temperature in the crosssection.

(5)Damping is not considered.

(6)Deflections and rotations are small before fluttering.

2.2.Basic equations

Bending and torsion of an open thin-walled beam are initiated mainly by the temperature gradients due to the absorbed heat flux.At the same time,the deformation also affects the incident angle of the heat flux.When the beam deforms,the absorbed solar heat flux is calculated as

where φ is the circumferential angle along the midline of beam cross-section;αsis the absorptivity of beam surface;S0is the magnitude of solar heat flux S0;α′∈ (0，2π)denotes the equivalent circumferential incident angle;ψ is the angle between S0and the deformed axis of the beam,and

θxis the torsion angle;θyiand θziare the bending angles of the centroid around y and z axis,respectively;δ is defined as

Based on Assumption(2),the beam temperature T( x，φ，t)is determined by

where c is specific heat;ρ is mass density;k is thermal conductivity;ε is the emissivity of beam surface;σ is the Stefan–Boltzmann constant.

Eq.(4)is a strong nonlinear equation,which is difficult to solve.However,it can be decomposed into two very simple equations by using the Fourier finite element method,16which approximates the temperature T( x，φ，t)as the sum of an average temperature Ta(x，t)and three perturbation temperatures:

Substituting Eq.(5)into Eq.(4)and integrating it over the cross-section with respect to φ,one can obtain two decoupled equations:

Eq.(6)is easy to solve because it is much simpler than Eq.(4).Upon obtaining the average temperature Ta(x，t)from Eq.(6),all perturbation temperatures can be solved by linear equation(7).

The nonuniform temperature distribution will result in thermal loads,which include two thermal bending moments MTyand MTzand a thermal bimoment BTas follows:

where E is the elastic modulus of beam;T0is the initial temperature of beam;αTis the thermal expansion coefficient;ω is the sectorial area.

According to Assumption(6),the change of thermal loads along beam length can be neglected.Thus,the vibration equations subject to these thermal loads are

with the boundary conditions:

where Dpis the torsional rigidity;Iyand Izare the moment of inertia around y and z axis,respectively;Iωis the sectorial quadratic moment;ρ is the density of the beam;v and w are the deflections of the shear center(yc,zc)in y and z direction,respectively;IAis the polar moment of inertia of the shear center;λ is defined asFor the open circular beam shown in Fig.1,

Let viand widenote the deflections of the centroid in y and z direction,respectively.They can be calculated by

According to Euler-Bernoulli beam theory,it yields

where χ = ?θx/?x.

Eq.(6),Eq.(7)and Eqs.(14)–(19)compose a set of coupled thermal-structural dynamic equations,which are difficult to solve due to these coupled terms and the nonlinearity of Eq.(6).However,the first Lyapunov method only investigates the linear approximations of these equations in their steady state for the stability analysis,so that it is not necessary to solve these equations directly.

2.3.Coupled steady state

As aforementioned,the solution of coupled steady state is generally required for the stability analysis and it cannot be obtained analytically.However,this coupled steady solution is slightly different from the uncoupled steady solution for a closed thin-walled beam,17in which the uncoupled solution is enough for stability analysis.Frustratingly,this is not true for an open thin-walled beam as illustrated by the numerical results shown in Section 3.Therefore,a simple iterative method has to be used to find an approximate solution.

Let t→ ∞ and x=l in Eqs.(6)–(19),and one can obtain the average temperature at iteration n(n=0,1,2,...,nmax):

Because Tachanges much more slowly than Tpi,19it can be regarded as a constant when Tpiis calculated.Therefore,

Consequently,

Then,one can update the following angles:

To start the iteration,the initial values can be set as

Through Eqs.(24)–(30),the approximate solutions of coupled steady state at the beam free end can be obtained within a small number of iterations.

2.4.Approximate solutions

Considering the first-order vibration mode of the beam,the deformations can be represented as

where N( x)and Φ( x)are the shape functions that satisfy the boundary conditions;V( t),W( t)and Θ( t)are the functions of t by using the method of separation of variables.

Substituting Eq.(31)into Eqs.(14)–(16)and noticing Eq.(20),one can obtain the following equations by using the Galerkin weighted residual method:

w h e r e

The above equations are related to the perturbation temperatures at x=l.As mentioned in Section 2.3,the average temperature can be regarded as a constant.Therefore,the perturbation temperatures can be solved from Eq.(7):According to Eqs.(33),(35)and(12),the deflection v is only related to the unknown variable Tp2(l，t),regardless of

Tp1(l，t)and Tp3(l，t).For the same reason,Tp2(l，t)is useless when we analyze the stabilities of deflection w and torsion.

Therefore,the stability of the deformations in different directions is discussed separately in the following.

2.5.Stability analysis

2.5.1.Sub-criterion A

To analyze the stability of the deflection v,one can rewrite Eqs.(33)and(36)in state space as

Eq.(38)is nonlinear because of the coupled term in matrix B.According to the first Lyapunov method,its asymptotical stability is determined by its linear approximation in steady state.Therefore,Eq.(38)is approximated about the steady state by using the first-order Taylor expansion as

According to Eqs.(2),(22)and(37),F can be calculated by

The characteristic polynomial of the matrix in Eq.(42)is

where s is the characteristic root;a0,a1and a2are corresponding coefficients.

A feasible shape function that satisfies Eq.(32)is

Substituting Eq.(47)into Eqs.(42)and(46),one obtains

According to the Routh–Hurwitz criterion,the stability conditions for a third-order linear system are

Obviously,a1＞0 and a2＞0 are always true.According to Eq.(48),a0＞0 is also satisfied because the term Izhk is much greater than FR4παTρch2l for most space beams.

Simplifying the last condition a1a2-a0＞0 in Eq.(51),one obtains

According to Eqs.(37)and(44),one can eventually get

Particularly,under the pure bending state in Y direction(around Z axis,α =0°or α =180°),θx= θy=0,and then this criterion becomes

2.5.2.Sub-criterion B

Compared to deflection v,the stability analysis of deflection in z direction and torsion is more difficult because these two deformations are coupled with each other,which leads to more unknown state variables.

The shape functions that satisfy Eq.(32)are

With these shape functions,Eq.(34)can be rewritten as

where

And according to Eqs.(11),(13)and(35),one obtains

Similar to Section 2.5.1,one can rewrite Eqs.(34)and(36)in the state space:

where the state variable vector is S=[W( t)

and

The linear approximation of Eq.(58)is

where

The parameters in Eq.(63)are

The characteristic polynomial of the matrix in Eq.(63)is

According to the Routh–Hurwitz criterion,the stability conditions for a sixth-order linear system are

where c1,d1,e1,f1and g1are the parameters of Routh table:

Eq.(72)is the stability condition for the z-direction bending and torsion coupling vibration.

2.5.3.Thermal- flutter criterion

The thermal- flutter criterion for this open thin-walled circular cantilever beam composes of the sub-criterion A given in Eq.(53)and the sub-criterion B given in Eq.(72).Unstable TIV will happen when any one of these two sub-criteria is violated.The sub-criterion A establishes the relationship between the incident angle of the solar heat flux and the stability of deflection v.Under the pure bending state in Y direction,this new criterion(Eq.(54))can degenerate into the existing criterion for a closed thin-walled beam17free of torsion and warping.The sub-criterion B establishes the relationship between the incident angle of the solar heat flux and the stability of deflection in Z direction and torsion.It is too complex to find a clear physical meaning as that of the sub-criterion A.However,it can be easily verified numerically.

3.Numerical results

In this section,numerical simulations based on the Fourier finite element method20will be conducted to obtain the dynamic responses of an open thin-walled circular cantilever beam subject to suddenly applied solar heat fluxes.The geometry dimensions and material properties are listed in Table 1.As shown in Fig.1,two incident angles are interested in:the normal angle θ0and the circumferential angle α.

Fig.2 depicts a typical vibration curve,from which one can easily identify that the vibration period is about 5 s.Therefore,the time step is set to 0.25 s in the following numerical simulations to ensure the numerical accuracy.

Figs.3(a)and(b)compare the uncoupled and coupled thermal-structural response of the bending angles θy, θzand the torsion angle θxat beam free end under solar fluxes of different incident angles α and θ0.It is clear that the uncoupled displacements are quite different from the coupled displacements,because the torsion angle θxhas great impact on the incident solar flux during the deformation of the beam.In order to get the coupled steady state values,the iterative method proposed in Section 2.3 with the maximum iterative number nmax=5 is utilized.The obtained results are θx=0.17 rad,θy=0.07 rad,θz=-0.03 rad for Fig.3(a)and θx=0.19 rad, θy=0.08 rad, θz=-0.01 rad for Fig.3(b),which are exactly the same as those from dynamic analysis.

As Section 2 emphasizes,when either sub-criterion A or sub-criterion B is violated,the thermal- flutter will happen.Based on this rule,one can plot the stable and unstable zones of the TIV of this beam in Fig.4.

Since this open thin-walled beam has ultra-low torsional stiffness,it is not strange that most cases are unstable in Fig.4.For example,when α =135°and θ0=0°,the steady--angles arerad.In this case,both sub-criterion A and sub-criterion B are violated.Accordingly,in the numerical simulation results depicted in Fig.5(a),all displacements are unstable.Thus,it verifies the prediction by the criterion.

Table 1 Geometry dimensions and material parameters.

Fig.2 Period of TIV when α=135° and θ0=0°.

Fig.3 Uncoupled and coupled thermal-structural responses at different incident angles.

Fig.4 Stable and unstable zones.

Fig.4 also implies a stable TIV when α =75°and θ0=30°.In this case,the steady angles are=0.21 rad,i=0.086 rad andi=-0.0016 rad,which satisfy both sub-criterion A and sub-criterion B.This conforms with the numerical simulation results depicted in Fig.5(b).

Fig.5 TIV at different incident angles.

An interesting case is pure bending state in Y direction.For example,when α=180°and θ0=20°,it is obvious that deflection w and torsion angle θxare equal to zero due to the symmetry of this problem.In this case,the sub-criterion A degenerates towhich is the same as the criterion for a closed thin-walled beam.17Since the steady anglesv could be stable according to the sub-criterion A.However,the sub-criterion B is violated,so that w and θxmust be unstable.These predictions are verified by the numerical results shown in Fig.5(c),in which all deformations are stable before 3000 s,then w and θxgradually diverge,and finally v is also unstable due to the influence of w and θx.This example demonstrates that the torsion has great impact on the stability of the TIV for open thin-walled beams.

4.Conclusions

This paper established a thermal- flutter criterion for an open thin-walled circular cantilever beam,which can be decomposed into two sub-criteria for the deflection v(sub-criterion A)and the coupled deformation of the bending state in z direction and torsion(sub-criterion B),respectively.In practice,thermal flutter happens when either the sub-criterion A or the subcriterion B is violated.The sub-criterion A can degenerate to the existing criterion for closed thin-walled beam free of torsion.However,the sub-criterion B shows that the torsion has great impact on the stability of the TIV for open thin walled beams that have ultra-low torsional stiffness.However,this criterion does not consider the structural damping,so that it gives a conservative prediction for practical structures.