Rigid–flexible–thermal coupling dynamics of a hub and multiplate system considering frictional contact

A geometric nonlinear modeling approach for strong rigid–flexible–thermal coupling dynamics of a hub and multiplate system considering frictional contact is proposed. Based on the absolute nodal coordinate formulation (ANCF), an ANCF thin-plate element with thermoelasticity is developed, where the temperature field is expressed with Taylor polynomials to yield heat-conduction equations. In contrast to the traditional coupling formulations, the influences of the attitude motion and structural deformation on the intensity of the solar radiation, the geometric nonlinearity of the plate as well as the frictional contact are taken into account. The frictional-contact formulations for a thin plate and a rigid body are presented, which can capture the stick–slip transition and address the multiple-point contact scenarios. To solve the strong rigid–flexible–thermal coupling equations, a novel numerical approach combining the generalized-α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\alpha $\end{document} method and the modified central-difference method is proposed. Two validations are performed to verify the proposed model, which proves the importance of considering the geometric nonlinearity and reveals the phenomena of thermally induced vibrations. Then, the thermal–dynamic coupling analysis for the satellite and solar-array multibody system in a thermal environment is carried out. The dynamic characteristics of the thermally induced vibration can be successfully revealed by the rigid–flexible–thermal coupling model. Furthermore, it is indicated that the influence of contact and thermal load on the nonlinear behavior of the solar-array deployment is essential, which demonstrates the feasibility of the proposed approach.


Introduction
Thermally induced vibrations of flexible spacecraft appendages may be initiated through the rapid changes of thermal loading during the orbital-eclipse transitions.The thermal gradient J. Liu liujy@sjtu.edu.cn in the thickness direction is caused by the difference of the heat flux applied on the upper and lower surfaces of the solar arrays, which may induce thermal deformation, thermal vibrations, and even influence the attitude motion of the satellite.Additionally, the thermaldynamic modeling of the spacecraft in many cases involves contact interaction between flexible appendages undergoing large deformations and overall motions.Therefore, an accurate rigid-flexible-thermal coupling model is essential to guide the practical engineering applications, and this study remains a worthy and challenging area.
Concerning the investigation of the thermal-structural coupling effect, Thornton and Kim [1] developed an analytical approach for determining the dynamic response of a flexible rolled-up solar array due to a sudden increase in external heating, in which the rotational motion of the spacecraft was not taken into account.Johnston and Thornton [2] presented an analytical model to investigate the effects of thermally induced vibrations of the flexible appendages on the rotational motion of the spacecraft, and then further carried out an experimental investigation on the thermal-structural coupling dynamics [3].Fan and Liu [4] developed a 2D rigid-flexible-thermal coupling model to investigate the thermally induced vibration of a hub-beam system.Shen et al. analyzed the thermal-structural coupling performance for a space thin-walled beam [5] and a spinning spacecraft with an axial boom [6], and further conducted a stability analysis for the thermally induced fluttering phenomena [7].Čepon et al. [8] proposed a new model for a coupled thermal-structural analysis of the bimetallic strip based on the ANCF.Liu et al. [9][10][11] carried out a thermal-structural analysis for a flexible spacecraft with a single or double solar panels, in which the coupling effect among attitude motion, structural deformation, and thermal loading are taken into account.Li et al. [12] investigated rigid-flexible-thermal coupling dynamics for a hub-beam system with a clearance joint considering torsional spring, latch mechanism, and attitude controller.
Since the solar arrays are thin plates, it is necessary to extend the thermal-structural coupling investigation to plate structures.Liu et al. [13] investigated the dynamic performance of a hub-plate system.The influence of the rotational motion on the intensity of the heat flux was included in the dynamic model.However, the temperature field expressed by the 3D solid element is complicated, and the geometric nonlinearity as well as the contact between the flexible bodies were not considered.
The Absolute Nodal Coordinate Formulation (ANCF) proposed by Shabana [14] has been widely used for simulation of a plate with large deformation.Since the transverse shear deformable of a thin plate can be neglected, based on Kirchhoff assumptions, Dmitrochenko et al. [15,16], Dufva [17], and Ren [18] parameterized the plate elements using slopes in the element midsurface direction only.Schwab et al. [19] and Sereshk and Mahmoud [20] compared the thin-plate elements against the plate element based on the conventional finite-element approach.Hyldahl et al. [21] also compared the convergence of rectangular thin-plate elements in different mesh configurations and load conditions.In ANCF, the gradient vectors are used to describe the nodal rotations instead of the rotation parameters, avoiding the singularity problem and the interpolation problem of the finite rotation variables [22][23][24].Hence, ANCF is suitable for the geometric nonlinear dynamics modeling of rigid-flexible coupled or rigid-flexible-thermal coupled multibody systems.Based on an ANCF shear-deformable plate element, Shen et al. [25] proposed a coupled thermalstructural model of a laminated composite plate.Cui and Yu proposed a novel method of the thermomechanical coupled analysis based on the unified description [26], and further developed a thermal integrated ANCF thin plate, which depicts the displacement and the temperature field integratedly [27].
Since the contact may occur during the deployment of a hub and multiplate system, it is necessary to consider the effect of frictional contact in the rigid-flexible-thermal coupling modeling.The key to modeling the contact dynamics of multibody systems lies in the treatment of contact-interface nonlinearities, which involves contact detection, contact discretization, enforcement of contact conditions, and a friction-force model [28][29][30][31].In recent years, some researchers have already developed several different formulations for modeling the frictional-contact problem in multibody systems.Konyukhov and Schweizerhof [32,33] developed a geometrically exact theory for various 3D solid-contact elements, including the normal, the tangential, and the rotational interactions.Yu et al. [34] established the frictionalcontact formulation between flexible and rigid bodies through the Hertz contact model and the velocity-based friction model.Shi et al. [35] presented a rotation-free shell formulation and an extended contact discretization for the dynamics of multibody systems with large deformations and frictionless contacts.Sun et al. [36,37] developed a new 2D segment-tosegment algorithm of contact dynamics based on an ANCF and mortar method, including both frictionless and friction cases.Gay Neto et al. [38] proposed a formulation to handle pointwise frictional-contact interaction of a beam-shell system using the degeneration of the local-contact problem.Tang and Liu [39] modeled the frictional contact in sliding joints with clearances, and generalized the regularized Coulomb friction law into an ANCF beam and rigid hole considering the stick-slip transition.Lei et al. [40] presented a strong coupling modeling method for the rigid bodies and the large deformed beams in contact interaction with the granular matter to simulate the flexible-brush sampling process.To avoid mutual penetration of the plate element with large deformation, a mixed method combining the node-to-surface, edge-to-surface, and surface-to-surface contact elements based on ANCF was proposed by the authors [41,42].Nevertheless, few studies have considered the frictional contact interactions in the rigid-flexible-thermal coupling multibody systems.
In this paper, a strong rigid-flexible-thermal coupling dynamic model for a hub and multiplate system is proposed considering the geometric nonlinearity and frictional contact.The remainder of the paper is organized as follows.In Sect.2, an ANCF thin-plate element with thermoelasticity is developed.Two frictional contact elements based on the ANCF thinplate element are presented in Sect.3. In Sect.4, the heat-conduction equations with solar radiation are derived, and a complete rigid-flexible-thermal coupling model is introduced.In Sect.5, a new numerical approach combining the generalized-α method and the modified central-difference method is proposed to solve the strong rigid-flexible-thermal coupling equations.Two validations are carried out to confirm the correctness of the proposed thermal-structural coupling model in Sect.6. Section 7 presents two interesting engineering numerical simulations of the satellite and solar-array multibody system in a thermal environment.Summary and conclusions are presented in Sect.8.

ANCF thin-plate element with thermoelasticity
A thin rectangular plate element with length a e , width b e , and height h, is shown in Fig. 1.As shown in the figure, O-XYZ and O e -X e Y e Z e represent the inertial coordinate system and the element coordinate system, respectively.
The global position vector of an arbitrary point of element i in the midsurface can be defined using the linear combination of the element-shape function matrix S and the global nodal coordinate vector q, which is given by r = S(x, y)q i , q i = B i q, where q i represents the element nodal coordinate vector, and x and y are the spatial coordinates defined in the element coordinate system, respectively.B i is a Boolean matrix.In the absolute nodal coordinate formulation, q i contains the global position vectors and slopes that are the spatial derivatives of the position vectors.The element nodal coordinate vector can be written as where In the classical theory of linear thermoelasticity, the thermal strain tensor is a linear relationship with temperature change [43,44].For a thin-plate structure, it is assumed that it causes only a positive strain but not a shear strain [45].Considering a plate element whose temperature is raised from the reference temperature T r at which strains and stresses are zero, to the temperature T , the thermal strain of the thin plate due to the temperature change T = T − T r is [43,46] given by where α x and α y are the coefficients of thermal expansion along the x and y directions, respectively.
Based on Kirchhoff assumptions, the virtual work done by the elastic force with thermoelasticity of the ANCF thin element i can be written as where D represents the matrix of modulus, and ε represents the Green-Lagrange strain vector of an arbitrary point of the plate, which is given by where ε and κ represent the inplane Green-Lagrange strain vector and the curvature vector, respectively, which can be given as where n is a vector normal to the midsurface of the plate and n is the magnitude of n.
For an arbitrary point of the plate, the temperature change T can be expressed with Taylor polynomials as T 0 (x, y, t) = T (x, y, 0, t), Considering that the thickness of the thin plate is small and the higher-order terms of Eq. ( 7) can be neglected, T can be reduced as Here, T 0 represents the temperature change of the midsurface of the thin plate and T 1 represents the temperature gradient along the thickness direction, which are obtained by solving the heat-conduction equations, as discussed later in Sect. 4.
Substituting Eqs. ( 5) and (9) into Eq.( 4), the virtual work done by the elastic force can be rewritten as Variation of ε 0 and κ leads to n = rx r y , Y = rx S y − ry S x , (14) in which rx and ry are the 3 × 3 skew-symmetric matrices corresponding to r x and r y , respectively.
Substituting Eq. ( 11) into Eq.( 10), the elastic force of the element is given by Accordingly, the elastic force can be written as Defining ρ as the mass density, g as gravitational acceleration, and μ as the structural damping coefficient, the variational equations of motion read where M = B T i ( ae 0 be 0 ρhS T Sdydx)B i is the constant generalized mass matrix, and C = μM is the damping matrix, and Q g = B T i ( ae 0 be 0 ρhS T gdydx) is the generalized gravitational force vector.

Frictional contact formulations for a thin plate and a rigid body
In this study, two forms of contact are taken into account, i.e., plate-to-plate contact and plate-to-rigid body contact.This section derives two frictional contact elements on the basis of the ANCF thin-plate element, which can address multiple-point contact scenarios with large deformations and overall motions.
The normal contact is formulated using the penalty method [28,30], and the tangential contact is modeled by the Coulomb friction law with a penalty regularization considering the stick-slip transition [29,32].The virtual work of the contact force is where g N is the normal penetration, is the contact area.u cf is the relative displacement vector of a contact pair, and f N , f T are the normal and tangential force vectors of a contact pair, respectively.Contact between the contact pair occurs only if g N ≤ 0.

Contact between plates
Using the "master-slave" contact-detection approach [29], the contact is checked between each integration points of the slave and the master surfaces, which is named as the STS contact element.As shown in Fig. 2, the nodal vector for the STS contact element can be given by The normal penetration can be defined as where h M and h S are the thicknesses of the master and slave surfaces, respectively.Following Eq. ( 14), the unit normal vector is obtained via the crossproduct of gradient vectors of the master surface, i.e., n = n/n.Here, the direction of the unit normal vector is perpendicular to the master surface and points to the slave node [42,47].Thus, the variation of the relative displacement vector of the contact pair SP is written as where S cf is the shape function matrix of the contact element, which can be defined as where S m and S s are the shape functions of the master and slave surfaces, respectively.ξ and η are the convective coordinates defined in the local surface-element coordinate system.Based on the penalty method [28,48], the normal contact force of the contact pair SP is formulated as where ε N is the normal penalty parameter.
In order to describe the stick-slip transition, it requires the history variables from the previous step to describe the incremental nodal displacement [39].Here, the points related to the previous step are represented by a subscript (•) a .As shown in Fig. 2, the contact point of the slave node S a on the master surface is P a in the previous step.Then, the incremental relative displacement vector of the contact pair SP is [42] where Accordingly, the incremental relative displacement vector on the current tangent plane is The return-mapping scheme [32] for the Coulomb friction law with a penalty regularization is applied to measure the friction force, it introduces the trial tangential force vector as where f T a is the real tangential force vector in the previous step, and ε T is the tangential penalty parameter.According to the yield function of the Coulomb friction law [29], the real tangential force of the contact pair SP can be determined as where μ s and μ d are the static and dynamic frictional coefficients, respectively.Then, the element residual vector of the STS contact element can be summarized as

Contact between a plate and a rigid body
As shown in Fig. 3, the contact is checked between each integration points of the slave-plate surface and the master rigid surface, which is named as the STRS contact element.Thus, the nodal vector for the STRS contact element can be given by The normal penetration can be obtained by where n = A c n is the unit normal vector of the STRS contact element, and A c is the direction cosine of the rigid body, and n is the unit normal vector of the rigid surface, which is defined in the body-fixed frame O c -X c Y c Z c .Then, the shape function matrix of the STRS contact element is where is the shape function of the slave rigid surface, ρ P is the position vector of the point P , which is defined in O c -X c Y c Z c , and K c is a function of the Cardan orientation angle vector of O c -X c Y c Z c with respect to the inertial frame [51].Similar to the STS contact element, the element residual vector of the STRS contact element can be obtained as Therefore, the generalized contact force vector related to Eq. ( 18) can be assembled as where B ST S e and B ST RS e are the connectivity matrices of the STS and STRS contact elements, respectively.n c 1 and n c 2 are the numbers corresponding to STS and STRS contact elements, respectively.

Heat-conduction equations with solar radiation
As shown in Fig. 4, the thermal loads due to solar radiation are applied on a thin plate with large deformation and overall motion.According to the first law of thermodynamics and the Fourier law of heat conduction [43][44][45], the variational heat-conduction equations for the three-dimensional continuum can be given by the authors' previous study [13], as follows where c is the specific heat of the material, k is the thermal conductivity, Q T is the intensity of the heat source, and q + , q − represent the heat flux of the upper and lower surfaces, respectively.
Since the distribution of solar-radiation input is affected by various factors, including its position, attitude, and configuration, the absorption of incident solar rays needs to be further evaluated by the heat flux [27].The solar unit vector is defined as s, the angle of incidence is α θ , and the initial angle of incidence is α 0 .Taking into account the effect of the shadow region, q + and q − can be given by where S 0 is the solar-radiation heat flux, c s is the solar coefficient related to the shadow region, α + m , α − m represent the effective solar absorptivity of the upper and lower surfaces, respectively, ξ + , ξ − represent the albedo view factors of the upper and lower surfaces, respectively, σ represents the Stefan-Boltzmann constant, ς + , ς − represent the emissivity of the upper and lower surfaces, respectively, and T + = T r + T + , T − = T r + T − represent the absolute temperature of the upper and lower surfaces, respectively.When the plate is in a shadow region, c s = 0.
Considering the influence of rigid-body motion and elastic deformation on the albedo view factors of the upper and lower surfaces, the albedo view factors are given by [2,13] where α + θ and α − θ are given by where the normal vectors are Since n+ and n− are closely related to the attitude motion of the rigid body as well as the elastic deformation of the plate, the present formulation considers the influence of the attitude motion and the structural deformation on the heat flux q + and q − , which is a complete rigid-flexible-thermal coupling model.
However, in the conventional formulation, with the neglect of the effects of rigid-body motion and elastic deformation on the heat flux, α + θ and α − θ are assumed to be constant, which can be expressed as In this case, it is not a complete rigid-flexible-thermal coupling model.
By means of Taylor polynomials, the three-dimensional temperature field of a thin plate can be described by the change of the midsurface temperature T 0 and the temperature gradient along the thickness direction T 1 .Using a finite-element method, T 0 and T 1 can be discretized as [49] T k = Np ik (k = 0, 1), (40) where p ik (k = 0, 1) represent the vectors of temperature variables of plate element i, and N represents the shape-function matrix.Defining p k (k = 0, 1) as the global vectors of temperature variables, the relation between the element vector and the global vector can be given by p ik = B ik p k (k = 0, 1).Substituting Eq. ( 40) into Eq.( 9), the absolute temperature in Eq. ( 35) can be expressed as Accordingly, the variational heat-conduction equations can be rewritten as where Considering that the temperature variable vector p is independent, the heat-conduction equations are given by 5 Dynamic equations for the rigid-flexible-thermal coupling multibody system T as the generalized coordinate vector of the rigid body, where r c is the position vector of the centroid defined in the inertial frame, and θ c is the Cardan orientation angle vector of O c -X c Y c Z c with respect to the inertial frame.In addition, the mass matrix M c and the generalized force vector Q c of the rigid body are detailed in [50,51].
Combining the constraint equations (q c , q, t) = 0, the rigid-flexible coupling dynamics equations [52] and the heat-conduction equations, the mixed differential-algebraic equations for the rigid-flexible-thermal coupling multibody system are given by where the generalized coordinate vector of the system is q d = [ q T c q T ] T , and the generalized mass matrix M d , damping matrix C d , and force matrix Q d take the form where n q is the number of the generalized coordinates of the plate, and λ is the vector of the Lagrange multiplier corresponding to the constraint equations.
A new numerical approach combining the generalized-α method and the modified central-difference method is proposed to solve the strong rigid-flexible-thermal coupling equations.The dynamic equations are solved by the generalized-α method [53,54], and the heat-conduction equations are solved by a modified central-difference method [49].Then, the Newton-Raphson algorithm is adopted to solve the combined nonlinear algebraic equations.The solution procedure can be divided into three steps, which are presented in detail in the Appendix.
Figure 5 shows the computational flowchart of the new numerical method proposed in this paper.Compared with the fourth-order Runge-Kutta method used uniformly in the authors' previous study [13], the computational efficiency of this proposed method is significantly improved.Additionally, compared with the traditional uncoupled method, this method ensures the accuracy of the full coupled simulation by solving the heat-conduction equations and the dynamic equations simultaneously.Considering the geometric nonlinearities (large deformation) and boundary nonlinearities (frictional contact), the variablestep-size control scheme [55] and the preconditioning strategy are employed to improve efficiency.

Validations
In this section, in order to verify the correctness of the proposed thermal-structural coupling model, simulation experiments of a plate clamped at two sides and a cantilevered plate applied with thermal loads are carried out, while the results from the commercial finiteelement software ABAQUS are taken as benchmarks.c = 921 J/(kg • K).In this simulation, the heat flux of the upper and lower surfaces are q + = q T cos α 0 and q − = 0, respectively, where q T = 1067 W/m 2 and α 0 = 0.
As shown in Fig. 7, the time histories of the midsurface temperature of point P and the temperature gradient T 1 = ∂T /∂z of point P are given, and the results are in good agreement with the comparative results provided by ABAQUS, which demonstrates that the temperature field can be calculated precisely by the proposed method.
To further investigate the geometric nonlinear effects, the deflections of point P in the z direction using linear and nonlinear models are compared in Fig. 8. Since the two edges of the plate are fixed, the thermal load leads to compressive thermal stress.It is found that for the nonlinear model, the thermal force causes the softening of the plate, which induces large transverse deformations.However, the results obtained by the linear model are quite different.Due to the neglect of the nonlinear stiffness matrices, the softening effect is not shown and the transverse deformation is small.Furthermore, the results of this paper coincide with the results achieved by ABAQUS, which validates the effectiveness of the proposed thermal integrated ANCF thin-plate element in terms of geometric nonlinearity.

Simulation of a cantilevered plate
As shown in Fig. 9, the dynamic simulation of a cantilevered plate applied with thermal loads is performed.The geometric properties and material data of the plate are the same as the numerical example 6.1.The time history of the temperature gradient of point P in the z direction is given in Fig. 10(a), which shows that the results are in accordance with the results obtained by ABAQUS.From the results of the deflections of point P in the z direction shown in Fig. 10(b), we can obtain that the results of the linear and nonlinear models coincide well, it can be concluded that for cantilevered plate without tip constraint, the use of a linear model is computationally accurate and efficient.The thermal-structural analysis for the cantilevered plate applied with thermal loads is further conducted through a comparative study of the thermal-structural coupling model and the constant heat flux model, where q T = 2525 W/m 2 and α 0 = 65 • , which ensures that q + = 1067 W/m 2 .The time histories of the temperature gradient and the deflection of point P in the z direction are shown in Fig. 11, respectively.As can be seen in the figures, due to the neglect of the elastic deformation of the constant heat-flux model, such a model cannot reveal the thermally induced vibration effect, therefore the vibration amplitude of the temperature gradient and the deflection is constant.On the contrary, since the heat flux varies with the deformation of the plate, the vibration amplitudes of the temperature gradient and the deflection obtained by the thermal-structural coupling model keep increasing.In summary, the proposed model takes into account the coupling effect of the heat flux and the elastic deformation, so that it is able to capture the phenomenon of thermally induced vibration, which is hardly achieved in the commercial finite element software ABAQUS.

Numerical simulations of a satellite and solar-array multibody system in a thermal environment
Thermal-dynamic analysis of the solar array is of vital importance to the safe operation of spacecraft.In this section, we present two numerical examples to demonstrate the effectiveness of the proposed rigid-flexible-thermal coupling model for the satellite and solararray multibody system in a thermal environment.The first numerical example focuses on comparing the different coupling models to reveal the dynamic characteristics of thermally induced vibration.The other numerical example is employed to analyze the deployment dynamic performance of the solar array considering contact and thermal loading.

Numerical simulation of a satellite and solar-array system in a thermal environment
A satellite and solar-array system applied with solar radiation is shown in Fig. 12.The flexible solar array is fully deployed, and its left edge is fixed to the rigid satellite.In this simulation example, the rigid satellite is a cube, and the solar array is an isotropic flexible  1. Defining the initial angel of incidence α 0 as the solar angle, the coordinate vector of the solar sunlight is given by s = sin α 0 0 − cos α 0 T .Initially, the satellite is in a static state, and the body-fixed frame of the satellite is parallel to the global coordinate system.The reference temperature is 273 K.

Dynamic performance of an onorbit satellite and solar-array system
To reveal the advantages of the rigid-flexible-thermal coupling model proposed in this paper for the satellite and solar-array system, two dynamic models are compared: the complete rigid-flexible-thermal coupling model (RFTC) and the rigid-flexible coupling model considering the thermal effect (RFCT).In the RFTC formulation, α + θ and α − θ are closely related to q t+ t d .However, in the RFCT formulation, α + θ and α − θ are assumed to be constant, which are equal to α 0 and π − α 0 , respectively, and thus the temperature results are not influenced by the rigid-body motion and the elastic deformation.
For a more intuitive description of the deformation, deflection is introduced.The deflection of an arbitrary point on the plate can be written as where r 0 is the origin of the body-fixed frame of the plate, and A p is the transformation matrix of the body-fixed frame with respect to the inertial frame.
The time histories of the deflection of the corner point P in the z c direction and the temperature gradient T 1 = ∂T /∂z of the corner point P are shown in Figs.13(a) and (b), respectively.As can be seen in the figures, due to ignoring the coupling effect of the solar heat flux, the rigid-body motion, and the elastic deformation, the RFCT model cannot reveal the thermally induced vibration.Accordingly, the vibration amplitude of the deflection is constant.On the contrary, the vibration amplitudes of the deflection and the temperature gradient obtained by the RFTC model keep increasing, which explains the phenomenon of the thermally induced vibration effect.
In order to further analyze the effect of thermal load on the motion of the satellite, Fig. 14 gives the time histories of the centroidal velocity of the satellite in the z direction and the angular velocity of the satellite in the y c direction.Due to the coupling of the translational and the rotational motion of the satellite and the elastic deformation of the solar array, the vibration amplitudes of the angular velocity and the centroidal velocity of the satellite obtained by the RFTC model also keep increasing.In summary, it reveals that the thermally induced vibration phenomenon is obvious in practical engineering and needs to be paid attention to.

Parameter analysis of thermally induced vibration
In this section, based on the proposed rigid-flexible-thermal coupling model, some characteristic parameters are investigated: the solar angle α 0 , the specific heat c of the plate, and the damping coefficient μ of plate on the thermally induced vibration effect.According to the change of the solar angle, the time history of the angular velocity of the satellite in the y c direction is shown in Fig. 15(a).It is indicated that with the decrease of the solar angle, the average of the deflection increases because of the increase of the intensity of the solar heat flux.However, the variation of the solar heat flux caused by the angular deformation of the plate and the rotation of the satellite becomes less significant, therefore, the thermally induced vibration due to the rigid-flexible-thermal coupling effect also decreases.
Figure 15(b) depicts the time history of the angular velocity of the satellite in the y c direction at different plate specific heats.As can be seen in the figure, due to the decrease of c and the increase of the thermal coefficient k/(ρc), the thermally induced vibration also weakens.
The time history of the angular velocity of the satellite in the y c direction under different plate damping coefficients is compared, as shown in Fig. 15(c).It can be seen that with the increase of the structural damping coefficient, the vibration amplitude of the deflection attenuates applied with the damping force, which indicates that the increase of the structural damping coefficient contributes to weakening the thermally induced vibration.

Dynamic performance of a maneuvering satellite and solar-array system
In this section, the nonlinear dynamic performance of the maneuvering satellite and solararray system applied with a thermal load is investigated.Additionally, the solar heat flux is applied on the plate, and the satellite is applied with a translational constraint, which is given by x c = 0, y c = 0, z c = 0, and a rotational constraint, which is given by ω x = 0, ω z = 0, and where t b = 200 s, t s = t b /8, ω s = π/50, which is shown in Fig. 16.
On the basis of the RFTC model proposed in this paper, the deformed configurations of the satellite and solar-array multibody system as well as the distribution of the von Mises stress of the solar array at some time steps are depicted in Fig. 17.Initially, the satellite rotates with angular velocity around the y c -axis, which in turn drives the motion of the solar array.In this situation, the solar array undergoes large deformations and overall motions.It is shown that the stress concentration mainly occurs at the roots of the solar array, and the stress wave gradually spreads toward the free edge of the solar array.
For the purpose of highlighting the advantages of the rigid-flexible-thermal coupling model for dynamic analysis of the maneuvering satellite and solar-array system with large deformations and overall motions, three dynamic models are compared: the RFTC model, the RFCT model, and the rigid-flexible coupling model without considering the thermal effect (RFC).In the RFC formulation, only rigid-flexible coupling effects are considered, while thermal effects are ignored.
Figure 18 gives the time histories of the deflection w of the corner point P in the z c direction and the temperature gradient T 1 = ∂T /∂z of the corner point P , respectively.As can be seen in the figure, in the RFC model, except for the acceleration and deceleration phases of ω y , the deflection of point P in the z c direction maintains a small-amplitude vibration due to the rigid-flexible coupling effect.For t < 200 s, the satellite undergoes periodic rotational motion in the y direction, which leads to the periodic deflection in the z c direction and the periodic variation of the temperature gradient in the thickness direction due to the complete rigid-flexible-thermal coupling effect, while the deflection and the temperature gradient obtained by the RFCT model shows neither periodic variation nor amplitude growing.Furthermore, it is shown that for t > 200 s, with the release of the applied angular velocity ω y , the amplitude of vibration obtained by the RFTC model keeps increasing slowly, which reveals the phenomenon of the thermally induced vibration.
To further analyze the thermal-dynamic performance of the maneuvering satellite and solar-array system, Fig. 19 gives the phase diagrams of point P in different models.From the phase diagram in the RFC model, it can be seen that the trajectory of the phase diagram exhibits a periodic motion without vibration under the action of ω y .In the RFCT model, the trajectory of the phase diagram can be viewed as the thermal vibrations superimposed on the basis of the RFC model.Since the thermal loads are assumed to be constant in the RFCT model, the thermal vibrations appear as low-frequency oscillations.However, in the RFTC model, it shows low-frequency periodic motion as well as the high-frequency thermally induced vibration, due to the consideration of the influences of rigid-body motion and elastic deformation on the intensity of solar radiation.Therefore, it is necessary to adopt the more precise rigid-flexible-thermal coupling model for dynamic analysis of maneuvering satellite and solar-array can capture richer nonlinear characteristics.

Deployment dynamics of a solar array in a thermal environment with frictional contact
The second example deals with the deployment of solar array in a thermal environment considering frictional contact.The simplified of the satellite and solar-array system in a thermal environment is illustrated by Fig. 20, which consists of four guide wires, eight solar panels, and eight rotational hinges.It is worth mentioning that the translational and attitude motion of the satellite is controlled during the deployment of the solar array, and it can be regarded as points B, D, B 1 , and D 1 are restricted to move on the guide wires ignoring the deformation of the guide wires.Since the direction of sunlight is along the negative z-axis, and the solar array are symmetrically distributed about the central satellite, it is sufficient to analyze the deployment of the solar panels on one side.O−XYZ is the inertial reference frame, which is consistent with the body-fixed frame of the satellite C − X c Y c Z c .The length of one side of the solar array is L P = 5 m, and the solar angle α 0 = 0, and for the other geometric and material properties refer to Table 1.The initial angle between where t 1 = 12 s, t 2 = 108 s, t 3 = 120 s and t e = 200 s.s 0 = 0.25L P cos(θ 0 /2) is the initial x coordinate, and a d = (L P − s 0 )/t 1 t 2 ensures that the deployment of the solar array ends with a fully deployed configuration, and the maximum velocity of the drive is defined as v m = v c = a d t 1 on the complete rigid-flexible-thermal coupling model proposed in this paper, the changes of configuration at some time steps are shown in Fig. 21, and the color maps represent the distribution of the von Mises The initial stress the folded solar panels is set to zero, and then the solar panels are gradually deployed applied with the driving constraints.Due to contact, it can be found that the stresses of panels 1 and 2 are greater than those of panels 3 and 4 at t = 30 s and t = 60 s.As the thermal loads continue, panels 3 and 4 experience greater thermal stress, which can be seen in the configurations at t = 90 s and t = 120 s.After the solar array is fully deployed, we can observe that the stress distribution of the four solar panels is relatively balanced, accompanied by obvious thermal deflection.
Figure 22 presents the time history of temperature gradients of the points P , B, and Q during the deployment of the solar array, and these test points can comprehensively reflect the heating conditions of the solar array.It can be seen that for t < 50 s, the temperature gradient at point P has a zero value, which is due to the fact that point P is in the shadow region at the time of contact and cannot receive solar radiation.Since one of the plates connecting the point B is in the shadow region at the initial stage, the temperature gradient rate at point B is lower than that at point Q.This why panels 1 and 2 experience greater thermal stress than panels 3 and 4 during the deployment of the solar array.Finally, the solar panels are fully deployed, and the temperature gradients of the three points are the same.Therefore, the influence of contact and thermal loads on the deployment process of the solar array is essential, which needs to be paid more attention to in engineering.

Influence of contact
In order to illustrate the necessity of contact, two models are employed.The first one considers the integration of the STS and STRS contact elements, and the other does not consider contact.For the sake of intuitiveness, the configuration of the satellite and solar-array system at t = 30 s is selected for analysis, as shown in Fig. 23.can be seen in Fig. 23(a) that the solar array can deploy properly without overpenetration in consideration of contact.On the contrary, it is clearly observed that the large penetration appears not only between the solar panels, but also between the satellite and the solar panels when the contact is not taken into account.In conclusion, for the deployment of the satellite and solar-array system, it requires a combination of STS and STRS contact elements to avoid mutual penetration, thereby ensuring computation accuracy.
With the view of further analyzing the effect of contact, Fig. 24 depicts the time history of coordinate x of points B and P obtained by the contact and no contact models.It can be seen that the values of the coordinate x of points B and P are less than 1, that is, the solar panels penetrate into the satellite, which violates the real situation.After the contact no longer occurs, the curves of the two contact models gradually tend to be consistent, and the solar array is fully deployed at t = 120 s.Therefore, we can conclude that the contact affects not only the motions of the system, but also the shadow regions, which is related to the intensity of the solar radiation.

Influence of thermal load
Since the thermal loads have a significant effect on the deployment of the solar array, two different thermal models are compared and analyzed in this section.Figure 25(a) presents the time history of coordinates z of point Q obtained by the thermal and no thermal models, respectively.By observing the coordinates z Q of the two models, it can be seen that when the solar panels 3 and 4 are fully deployed, the thermal model excites high-frequency oscillations with gradually decreasing amplitudes.Combined with the configuration at t = 55 s in Fig. 25(b), it can be seen that since point B of the no thermal model moves later than that of the thermal model, point Q has a greater vertical velocity as it passes the horizontal position, resulting in a greater offset magnitude.In addition, the final configurations of the thermal and no thermal models at t = 200 s are also presented in Fig. 25(b), in which we can observe that obvious deflections of the thermal model occur at points P and Q due to the effects of the boundary constraints and the thermal loads, which cannot be obtained by the no thermal model.To further analyze the thermal expansion and deflection, the time history of the deflection of point P in different models is given in Fig. 26.Here, u, v, and w denote the deflection of point P in the x, y, and z directions, respectively.The values of u and v remain zero in the no thermal model, but eventually exhibit positive expansion in the thermal model.By comparing the curve of w in the two models, it can be observed that when the solar array is fully deployed at t = 120 s, the vibration amplitude of the thermal model is greater than that of the no thermal model, and the decay time is longer.From the results in the figure, it is necessary to pay attention to the phenomenon of thermal deflection during the deployment of the solar array.

Comparison of different drive speeds
In engineering, the effect of drive speed on the satellite and solar-array system is another issue of concern.In this section, we focus on comparing the results of three different drive speeds acting on the solar array in a thermal environment.As illustrated in Fig. 20, the points D and D 1 are subjected to a prescribed translational displacement along the x direction.As shown in Fig. 27, the drive-speed curves in different cases have the same displacement at t = t e , in which the black curve represents the basic model, and its expression refers to Eq. (55). Figure 28 plots the dynamic response of point P at three different drive speeds.It can be seen in Fig. 28(a) that the deployment times of the solar array in the three cases are significantly different, and the changes of x P are relatively stable.From the results in Fig. 28(b), we can observe that there are two obvious high-frequency vibrations during the deployment of the solar array.The first occurs when the solar panels 3 and 4 are fully deployed, another occurs when the solar panels 1 and 2 are fully deployed.Since the thermal effect is positively related to the action time of the thermal load, the thermal stress increases with decreasing drive speed, resulting in a greater amplitude of oscillation, especially when the plates are fully deployed.

Conclusions
A strong rigid-flexible-thermal coupling dynamic model for a hub and multiplate multibody system with large deformations and frictional contact interactions is proposed.An ANCF thin-plate element with thermoelasticity is developed, and the temperature field is expressed with Taylor polynomials to gain the heat conduction equations.Accounting for the influences of the attitude motion and structural deformation on the intensity of the solar radiation, a complete rigid-flexible-thermal coupling model is created.Subsequently, the STS and STRS frictional contact elements are developed, which can effectively avoid the mutual overpenetration and capture the stick-slip behavior.Furthermore, a new numerical approach combining the generalized-α method and a modified central-difference method is proposed to solve the strong rigid-flexible-thermal coupling equations, which improves the accuracy of the full coupled simulation by solving the heat-conduction equations and the dynamic equations simultaneously.
By setting comparative tests in ABAQUS, the correctness of the proposed thermalstructural coupling model is validated, and the necessary of the geometric nonlinearity is proved.Moreover, the advantages of the proposed model in capturing thermally induced vibration are revealed.To further demonstrate the effectiveness of the proposed rigidflexible-thermal coupling method, two numerical simulations of the satellite and solar-array multibody system in a thermal environment are conducted.The first example compares the results of three dynamic models: RFTC model, RFCT model, and RFC model.Due to the consideration of the elastic deformation as well as the rigid-body motion on the intensity of heat flux, RFTC model can reveal the phenomena of the thermally induced vibrations, which cannot be achieved by the other two models.Parameter analysis illustrates that the thermally induced vibration effect can be reduced through the decrease of the solar angle and the specific heat or the increase of the damping coefficient.Then, the nonlinear coupling behavior of the maneuvering satellite and solar-array system can be revealed by the proposed RFTC model.The second example is employed to analyze the deployment dynamic performance of the solar array in a thermal environment with frictional contact.It is concluded that the contact affects not only the motions of the system, but also the intensity of the solar radiation because of the shadow regions.Due to the effects of the boundary constraints and the thermal loads, obvious thermal deflection and high-frequency vibration can be induced, especially when the plates are fully deployed.Since the thermal effect is positively related to the action time of the thermal load, the amplitude of thermal vibration increases with the decreasing of the drive speed.

Appendix. Solution scheme for the rigid-flexible-thermal coupling equations
The solution scheme for the rigid-flexible-thermal coupling equations is divided into the following three steps:

• Step1: Prediction of the initial value
Assuming that qt+ t d = 0, we combine Taylor's formula, the generalized-α method, and a modified central-difference method to predict a t+ t d , q t+ t d , qt+ t d , λ t+ t , p t+ t , as follows: where and the numerical parameters ρ ∞ , θ ∈ [0, 1].The numerical parameters in this paper are taken as ρ ∞ = 0.6, θ = 0.5.

• Step2: Transformation into nonlinear algebraic equations
The dynamic equations can be transformed into nonlinear algebraic equations using the generalized-α method [ The heat-conduction equations in Eq. ( 48) are transformed into nonlinear algebraic equations by a modified central-difference method [49] Kp t+ t − Q(q t+ t d , p t+ t ) = 0, where • Step3: Solution of the combined nonlinear algebraic equations According to Eqs. (A8) and (A11), the combined nonlinear algebraic equations are rewritten as The Newton-Raphson algorithm is used to find a solution of Eq. (A14).At an iteration step k, the following equation is solved for a correction q t+ t (k) : An improved estimate can be obtained as q t+ t (k+1) = q t+ t (k) + q t+ t (k) , k = 0, 1, . . ., (74) until the precision condition q t+ t (k) ≤ δ is satisfied, in which δ is the precision error.

Fig. 1 A
Fig. 1 A thin rectangular plate element of ANCF (Color figure online)

Fig. 2
Fig. 2 Geometry and kinematics of the STS contact element (Color figure online)

Fig. 3
Fig. 3 Geometry and kinematics of STRS contact element (Color figure online)

Fig. 4 A
Fig. 4 A thin plate with large deformation and overall motion under solar radiation (Color figure online)

Fig. 6 AFig. 7
Fig. 6 A plate clamped at two sides in a thermal environment (Color figure online)

Fig. 8 Fig. 9 AFig. 10
Fig. 8 Time history of the deflection of point P in the z direction (Color figure online)

Fig. 11
Fig. 11 Time histories of (a) temperature gradient; (b) deflection of point P in the z direction (Color figure online)

Fig. 12 ATable 1
Fig. 12 A satellite and solar-array system (Color figure online)

Fig. 13
Fig. 13 Time histories of (a) deflection of point P in the z c direction; (b) temperature gradient of point P in the z direction (Color figure online)

Fig. 14 Fig. 15
Fig. 14 Time histories of (a) centroidal velocity of the satellite in the z direction; (b) angular velocity of the satellite in the y c direction (Color figure online)

Fig. 16
Fig. 16 Time history of ω y (Color figure online)

Fig. 17 Fig. 18
Fig. 17 Motion of the satellite and solar-array system at some time steps (Color figure online)

Fig. 21
Fig. 21 Configurations and the von Mises stress of the satellite and solar-array system at some time steps (Color figure online)

Fig. 22 Fig. 23
Fig. 22 Time history of temperature gradients of points P , B, and Q (see Fig. 20) in the z direction (Color figure online)

Fig. 24
Fig. 24 Time history of coordinates x of points B and P in different contact models (Color figure online)

Fig. 26 Fig. 27
Fig. 26 Time history of deflection of point P in different thermal models (Color figure online)

Fig. 28
Fig. 28 Time histories of (a) displacement; (b) velocity of point P in the x direction in different cases (Color figure online) 53]