[BACK]
images Computer Modeling in Engineering & Sciences images

DOI: 10.32604/cmes.2021.016475

ARTICLE

Thermally Induced Vibration Analysis of Flexible Beams Based on Isogeometric Analysis

Jianchen Wu1, Yujie Guo1,* and Fangli Wang1,2

1College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing, 210016, China
2School of Mechanical Engineering, Jinling Institute of Technology, Nanjing, 211169, China
*Corresponding Author: Yujie Guo. Email: yujieguo@nuaa.edu.cn
Received: 08 March 2021; Accepted: 26 May 2021

Abstract: Spacecraft flexible appendages may experience thermally induced vibrations (TIV) under sudden heating loads, which in consequence will be unable to complete their intended missions. Isogeometric analysis (IGA) utilizes, in an isoparametric concept, the same high order and high continuity non-uniform rational B-splines (NURBS) to represent both the geometry and the physical field of the structure. Compared to the traditional Lagrange polynomial based finite element method where only C0-continuity across elements can be achieved, IGA is geometrically exact and naturally fulfills the C1-continuity requirement of Euler–Bernoulli (EB) beam elements, therefore, does not need extra rotational degrees-of-freedom. In this paper, we present a thermally induced vibration analysis framework based on the isogeometric method where thermal and structural behaviors are coupled. We fully exploited the higher order, higher continuous and geometric exactness of the NURBS basis with both benchmarks and sophisticated problems. In particular, we studied the thermally induced vibrations of the Hubble Space Telescope (HST) solar panel where main factors influencing thermal flutters are studied, and where possible improvements of the analytical reference methods are discussed. Additionally, thermally induced vibrations of the thin-walled lenticular tubes are studied and two new configurations of the tube are proposed to effectively suppress the thermally induced vibrations. Numerical examples of both benchmarks and sophisticated problems confirm the accuracy and efficiency of the isogeometric analysis framework for thermally induced vibration analysis of space structures.

Keywords: Thermally induced vibration; thermal flutter; radiation heat transfer; isogeometric analysis; thermal structural coupling

1  Introduction

Spacecraft in orbit periodically moves in and out of earth’s dark and sunny regions. The temperature field of the spacecraft changes rapidly at the day-night or night-day transitions due to the combined effects of the solar heat flux, the earth-emitted heat flux and the earth albedo heat flux [1]. Flexible appendages, commonly used in spacecraft, are known to have low stiffness and low frequencies, therefore, they are prone to vibrate even under small excitation forces. The typical excitation load in the flexible appendages is the time-dependent bending moments generated by the temperature gradients within the structure. Since the temperature field and the displacement field are coupled, thermally induced vibrations may happen. This phenomenon needs special attention especially for the flexible space structures, since unstable oscillation or thermal flutter may occur under certain conditions [2] which will degrade the flight quality of the spacecraft and affect its normal operations. A famous example in history is the Hubble Space Telescope (HST) launched in 1990 and it was suffered from the spherical aberration and a pointing “jitter” due to the thermally induced vibrations [3,4]. The Hubble Space Telescope solar arrays were then replaced by the US space shuttle “Endeavour” in 1993.

Analytical study of the thermally induced vibrations of a simply supported rectangular beam was first investigated by Boley in 1956 [5], and then extended to the structures with more complex shapes, such as plates and shells [68], etc. In [5], a dimensionless parameter B was proposed to determine the maximum deflection of the beam. This parameter is directly proportional to the beam’s natural frequency and the characteristic thermal response time, hence simplifies the analytical expressions. In addition, it is worth to mention that, in Thornton et al. [8] developed an analytical approach for determining the thermal-structural response of a flexible solar array and established the stability criterion to judge whether thermal flutter occurs. For experimental studies, torsional thermal flutter of booms with open cross section was observed and compared to the theoretical results in [9]. Then, Rimrott et al. [10] successfully triggered torsional thermal flutter of beams in experiments. Recently, Su et al. [11] conducted an experimental study of the thermally induced vibration of space boom structures and demonstrated the validity of theoretical and numerical methods to analyze such problems.

Despite the various analytical methods proposed in the literature [58,12], the practical space structures are often too complex to cope with for the analytical methods. Therefore, numerical methods such as finite elements are developed to study the thermal-structural coupled behaviors of the space structures. Manson [13] proposed, for the first time, a finite element model for the thermally induced vibration analysis of beams and plates. Later on, more complex situations, such as temperature-dependent material properties [14], composite and smart materials [15,16], design optimizations [17,18] and large-scale space structures [19,20] are considered in the finite element model. Recently, Liu et al. [21] developed a rigid-flexible-thermal coupling model to study the coupling effect among attitude motion, structural behaviors and the thermal loading of the spacecraft. Regarding different numerical methods, Shen et al. [22] developed an ANCF (Absolute Nodal Coordinate Formulation) method to analyse the thermally induced vibrations of flexible beams including large rotations. The ANCF method was further extended to study the thermal shock induced dynamics of a deploying boom in [23]. We note that, for thermally induced vibration analysis of beam structures, most of the researches are focusing on the circular cross sections, more complex shapes such as lenticular shapes are rarely considered.

Isogeometric analysis (IGA) [24] is an emerging computational approach aiming at bridging the gap between computer-aided design (CAD) and computer-aided engineering (CAE). Compared to the classical finite element method which uses C0-continuous Lagrange polynomials as basis functions, IGA employs the non-uniform rational B-spline (NURBS) basis functions to represent both the geometry and physical fields which allows for an exact geometric description and provides higher order and higher continuity basis functions for analysis. Due to its superior properties, IGA has been applied to a wide range of areas [2532], especially in the field of thin-walled structures, e.g., in the modeling of Euler–Bernoulli beam elements [3335] and Kirchhoff–Love shell elements [3641] where at least C1-continuity of the basis are needed. For isogeometric thermal analysis, the steady state and transient heat transfer analysis of solids are investigated in [42,43] where superior convergence rates of IGA are achieved. For isogeometric structural dynamic analysis, different aspects have been investigated, e.g., the free vibration problems [4446] and the nonlinear dynamics of structures [34,47], etc.

In this paper, thermally induced vibrations of beam-like structures in the framework of isogeometric analysis are investigated for the first time, where the higher order, higher continuous and geometrically exact properties of the isogeometric analysis are fully exploited. The coupling between thermal and structural responses are implemented to capture thermal flutter behaviors of the beam-like structures. In particular, a Hubble Space Telescope solar array assembled from three main parts are studied and the main factors influencing thermal flutters are investigated in detail. We also compared IGA results with the analytical reference solutions and discussed the possible improvements of the analytical methods. Additionally, we studied the thermally induced vibration responses of thin-walled lenticular tubes with different cross section shapes, which are widely used as supporting frames of solar sail due to its highly accurate and repeatable deployment properties [48]. Based on these studies, we propose two new configurations of the lenticular tube, where thermally induced vibrations can be significantly suppressed. We note that, with isogeometric analysis, various cross section shapes can be modeled exactly.

The paper is organized as follows: in Section 2 we provide a brief summary of the NURBS basis functions and NURBS geometries. In Section 3 we describe in detail the basic equations governing the thermal and structural behaviors of beam structures and the corresponding isogeometric discretization. In Section 4, we test our method with several numerical examples to reveal the method’s capabilities. Finally, we summarize the main aspects and findings and draw conclusions in Section 5.

2  NURBS Basis and Geometries

In this section, we first summarize the basic properties of the B-spline and NURBS basis in a brief manner, we then introduce the constructions of the NURBS geometries and its derivatives which are frequently used in this paper.

A NURBS curve of order p is defined as [49]:

C(u)=i=1nRi,p(u)Bi (1)

where Ri,p(u) are the NURBS basis functions with the parametric coordinates u, where Bi=(xi,yi,zi) are the control points of the NURBS curve in Cartesian coordinates, and where n is the total number of control points. The NURBS basis function Ri,p(u) can be represented as:

Ri,p(u)=Ni,p(u)ωii^=1nNi^,p(u)ωi^ (2)

where ωi is the weights of the i-th control points, and where Ni,p(u) is the i-th B-spline basis function of order p. The B-spline bases are defined through a non-decreasing knot vector:

Ξ={u1,u2,,un+p+1},with uiui+1(i=1,2,,n+p) (3)

where ui(i=1,,n+p+1) is the i-th knot value. The knot vector Ξ is said to be open when the first and the last knots are repeating p+1 times.

Based on a given knot vector Ξ, the B-spline basis functions can be obtained through the following Cox-de-Boor recursion formula [50]:

Ni,0(u)={1 if uiu<ui+10 otherwise (4)

Ni,p(u)=uuiui+puiNi,p1(u)+ui+p+1uui+p+1ui+1Ni+1,p1(u) (5)

Fig. 1 shows an example of cubic B-spline bases with the knot vector defined as Ξ={0,0,0,0,1,2,2,3,3,3,4,4,4,4}. It can be found that, for open knot vectors, the B-spline basis functions are interpolatory at both ends of the knot vector, while at the internal knots with multiplicity of k, the B-spline basis are Cpk-continuous. For example, in Fig. 1, the B-spline bases are C2, C1 and C0 continuous at the internal knots u=1,2,3, respectively. This unique property allows for a direct operation on the inter-element continuity of the NURBS basis, while for the traditional Lagrange polynomials, only C0-continuity across elements can be obtained.

images

Figure 1: 1D cubic B-spline shape functions Ni,3(u)=(i=1,,10) across an open knot vector of four knot-span-elements

A NURBS surface can be constructed based on the tensor product of one-dimensional geometries as:

S(u,v)=i=1nj=1mRijpq(u,v)Bij (6)

where Rijpq(u,v) are the two-dimensional NURBS basis of order p and q in the u and v parametric directions, respectively. They are defined as:

Rijpq(u,v)=Ni,p(u)Mj,q(v)ωiji^=1nj^=1mNi^,p(u)Mj^,q(v)ωi^j^ (7)

where Ni,p(u) and Mj,q(v) are the one-dimensional B-spline basis, where ωij are the weights associated to the (i,j)-th NURBS basis.

The derivatives of a NURBS curve can be simply obtained by the following relations:

C(k)(u)=i=1nRi,p(k)(u)Bi (8)

where the superscript (k) denote the k-th derivatives w.r.t. the parametric coordinate u. Based on Eq. (8), one can easily get the tangent vector of a NURBS curve by taking its first derivatives C(1)(u). The above rule also applies to the NURBS surfaces and solids.

3  Isogeometric Thermally Induced Vibration Analysis of Beam Structures

In this section, we first present the fundamental equations governing the thermal behaviors of beams and its isogeometric discretization in 3.1. We then introduce the basic equations for structural analysis of the beam structures based on the isogeometric analysis in 3.2. Lastly, a coupled scheme between thermal and structural responses are presented in 3.3.

3.1 Isogeometric Thermal Analysis of Beam Structures

We consider a thin-walled tube exposed to the solar heat flux S0, as shown in Fig. 2, where φ is the solar incident angle, where S=S0cosφ represents solar heat flux perpendicular to the tube axis. Various cross-sectional shapes can be assigned to the tube, e.g., circular, rectangular or lenticular shapes. Here h is the thickness of the tube. We assume that the tube is very thin such that the temperature gradient along the wall thickness can be ignored. Other adopted assumptions are: no heat conduction along the length of the tube; the heat convection is neglected due to the high vacuum space environment; the radiation heat transfer on the inner wall of the tube is not included.

images

Figure 2: Thin-walled tube subjected to solar heat flux with different cross section shapes

images

Figure 3: Heat conduction in the cross section of the tube: (a) problem set up, (b) line of incident for a general B-spline curve

Since the heat conduction along the beam axis is neglected, the problem considered here is simplified to a 2D plane problem where only the cross section is studied. Consider a differential element ds in the cross section of the tube under thermal equilibrium, cf. Fig. 3a, where s is the arc-length coordinate along the line of centroids of the cross-section. Here we denote the curve connected by the points located at the middle of the tube’s wall as line of centroids, which is in contrast to its traditional definition along the beam axis. In Fig. 3a, qs is the absorbed heat flux, qσ is the thermal radiation from the outer wall of the tube to space, and qin and qout are the heat fluxes conducted into and out of the tube element, respectively. The heat balance equation of the differential element can be formulated as:

Ttksρc2Ts2+σερchT4αsρchqs=0 (9)

where ks is the thermal conductivity along the arc length direction, where ε and αs are the emissivity and absorptivity of the tube’s surface, respectively, where σ is Stefan–Boltzmann constant, c is the specific heat, where ρ is the density of the tube. As we mentioned before, Eq. (9) can be applied to arbitrary cross sections of thin-walled tube. For annular shapes, cf. Fig. 2, the absorbed heat flux is expressed as qs=Sδ(θ)cos(θ), where δ(θ) is a function of θ determining which side of the tube is under radiation and can be written as:

δ(θ)={1,π/2<θ<π/20,π/2<θ<3π/2 (10)

For general cross-sectional shapes, where the line of centroids of the cross section is represented by a B-spline curve C(u), cf. Fig. 3b, the absorbed heat flux qs is written as:

qs=Sδ(θ)cos(θ1)=Sδ(θ)sin(θ2)=Sδ(θ)sin<S,C(1)(u)> (11)

where C(1)(u) is the first derivative of the B-spline curve w.r.t. the parametric coordinate u and can be obtained from Eq. (8).

Applying the method of weighted residuals (MWR) and integrating Eq. (9) by parts, we obtained a weak form of the governing equation:

AρcTtT˜dA+AksTsT˜sdA=A(αshqsσεhT4)T˜dA (12)

where A is the cross-sectional area of the thin-walled tube, where T˜ is the test function, and where we assume the tube is under adiabatic boundary conditions.

Following the concept of isogeometric analysis, the geometry and the temperature field are discretized with the same NURBS basis functions:

T=T(s,t)=i=1nNi,pTi=NTT (13)

T˜=T˜(s,t)=i=1nNi,pT˜i=NTT˜ (14)

where Ti denote the temperature unknowns of control points, where N and T are the vectors collecting all the basis functions and unknowns, respectively.

Substitute Eqs. (13)(14) into (12) and after math manipulations, we obtain the transient heat transfer equation:

CdTdt+KT=Rq+Rσ (15)

where K and C are the conductance and capacitance matrices and can be written, respectively, as:

K=Aks(Ns)TNsdA (16)

C=AρcNTNdA (17)

where N/s can be further expressed as:

Ns=Nuus (18)

In Eq. (15), Rq and Rσ are the load vectors arising from specified surface heating and surface radiation, respectively, and can be written as:

Rq=AαshqsNTdA (19)

Rσ=AσεhT4NTdA (20)

Eq. (15) is a nonlinear transient heat conduction problem due to the nonlinear term Eq. (20). Here we combine the implicit one-parameter θ˜ time integration scheme and the Newton-Raphson method to solve Eq. (15) iteratively. Here we use θ˜ instead of θ to represent the time integration parameter. The time integration scheme can be formulated as [51]:

K¯Tt+Δt=R¯t+Δt (21)

where the subscript denotes the time step, where Δt is the time step increment, and where K¯ and R¯t+Δt are formulated, respectively, as:

K¯=θ˜K+1ΔtC (22)

R¯t+Δt=[(1θ˜)K+1ΔtC]Tt+(1θ˜)Rt+θ˜Rt+Δt (23)

where Rt is the load vector at time t:

Rt=(Rq+Rσ)|t (24)

During each time step, the Newton-Raphson method was adopted to solve the nonlinear equations [49]

JmΔTm+1=Fm (25)

Tm+1=Tm+ΔTm+1 (26)

where the superscript m is the index of iteration, where J is the Jacobian matrix and F is a residual load vector which are represented as:

Jm=K¯+ΔRm (27)

Fm=K¯TmR¯m (28)

where ΔR is the contribution from the temperature-dependent heat load vector, as:

ΔRm=4Aσεh(Tm)3NTNdA (29)

Eqs. (25) and (26) are computed iteratively until the residual F reaches the predefined tolerance.

3.2 Isogeometric Structural Analysis of Beams

The Euler–Bernoulli beam model is used in this paper to study the structural behaviors of the thin-walled tube. The governing equation of the beam is written as:

EI4wz4+2MTz2+ρA2wt2=0 (30)

where w is the beam’s deflection, where E and I are the modulus of elasticity and the area moment of inertia, respectively, and where MT is the thermal moment, defined as:

MT=AEαTΔTydA (31)

where αT is the coefficient of thermal expansion, where ΔT is the temperature gradient of the beam’s cross section which can be obtained as:

ΔT=T(s,t)Tave(t)=T(s,t)0LsT(s,t)dsLs (32)

where Tave(t) is the average temperature of the cross section at time t, and where Ls is the length of the line of centroids of the cross section. Taking an annular cross section for example, the inner and outer radius of the cross section are R1 and R2, respectively, then the length of line of centroids is π(R1+R2).

Applying the method of weighted residuals (MWR) and integrating Eq. (30) by parts twice yields the weak form of the governing equation:

0LEI2wz22w˜z2dz+0LMT2w˜z2dz+0LρA2wt2w˜dz=0 (33)

where L is the length of the thin-walled tube, where w˜ is the test function of the displacement field.

Similar to the isogeometric thermal analysis, the displacement field of the beam is discretized with the same NURBS basis as the geometry:

w=w(z,t)=i=1nNi,pwi=NTw (34)

w˜=w˜(z,t)=i=1nNi,pw˜i=NTw˜ (35)

where wi is the displacement unknowns at the i-th control point, and where w is the vector collecting all the displacement unknowns.

Substitute Eqs. (34), (35) into (33) results in the discretized form of governing Eq. (33):

Md2wdt2+Ksw=FT (36)

where Ks and M are the stiffness matrix and consistent mass matrix, respectively, they are defined as:

Ks=0LEI(d2Ndz2)Td2Ndz2dz (37)

M=0LρANTNdz (38)

where FT in Eq. (36) is the equivalent nodal load vector induced by the thermal moment MT:

FT=0LMT(d2Ndz2)Tdz (39)

Eq. (36) can be solved iteratively by the Newmark algorithm [51]. In the following numerical examples, the parameters of the Newmark algorithm are selected as: α=0.25, β=0.5, which ensure an unconditionally stable algorithm.

3.3 Thermal Structural Coupling

The temperature and displacement fields of the beam are coupled with each other. The beam will deform under thermal moments while, at the same time, this deformation will change the incident angle of the sun which has direct influence on the temperature field of the beam. Here we modified the heat load vector Rq in Eq. (19) to take into account the coupling effect:

Rq=AαsS0cos(φw/z)hδ(θ)cos(θ1)NTdA (40)

where w/z is the slope of the beam. In isogeometric analysis, w/z is expressed as:

wz=(NTw)z=NTzw=NTuuzw (41)

Fig. 4 illustrates the change of incident angle induced by the thermal-mechanical coupling. In this paper, the sequential coupling method are used to compute the thermally induced vibrations of beam structures where, in a single time step, we first compute the temperature field of the beam, we then apply the obtained thermal moment for vibration analysis. In all the numerical examples, we adopt a time step equals to 0.1 s.

images

Figure 4: The change of solar incident angle in coupled thermal-structural analysis

4  Numerical Examples

In this section, we present a number of numerical examples to demonstrate the reliability, accuracy and robustness of the proposed method. We start with a simply supported rectangular beam benchmark example, cf. 4.1 to verify the accuracy of the proposed method. We then consider a more realistic example, a solar array of Hubble Space Telescope (HST) to study the influential factors of the thermal flutter phenomenon in 4.2. We compare our IGA results with the reference solutions which is less accurate due to the approximations introduced in the reference model. In the last example, we studied the influence of the cross-section shapes on the thermally induced vibrations of the thin-walled tube. Particularly, we proposed a new type of thin-walled lenticular tube in 4.3, where thermally induced vibrations can be effectively suppressed.

4.1 Thermal Vibration of a Simply Supported Rectangular Beam

In this example, a simply supported rectangular beam subject to a constant heat flux Q on the top surface (y=h/2) is studied. The geometry descriptions and boundary conditions of the beam are shown in Fig. 5. The bottom surface (y=h/2) of the beam is adiabatic. We adopt two separate discretization of the beam for thermal and structural analysis to take into account the difference between the two governing equations, namely the first order and the second order Partial Differential Equations (PDEs). For thermal analysis, the beam is discretized with 4 elements with order of basis p=2, while for structural analysis, 8 elements with order of basis p=3 is used.

images

Figure 5: Simply supported rectangular beam under surface heat flux

The analytical solution of the beam is given by Boley [5] as:

V=[(ξ2ξ)/2]mT+n=1,3,5sinnπξn3π3[π28B2n2sinn2π2B2τj=1,3,5ej2π2τ+(j/nB)2sinn2π2B2τcosn2π2B2τj4+n4B4] (42)

where ξ[0,1] is the normalized coordinates in the length direction, where mT is the non-dimensional thermal moment, which is defined as:

mT=π4kMT192EIQαT (43)

where k is the thermal conductivity of the beam.

The non-dimensional time τ and parameter B in Eq. (42) are defined, respectively, as:

τ=κth2 (44)

B=hLκ(EIρA)14 (45)

where κ=k/(ρc) is the thermal diffusivity. The parameter B defined in Eq. (45) is non-dimensional and represents the thermal and mechanical properties of the material. In this example, we set B=1.

Fig. 6 shows a comparison between the mid-span displacements of the IGA model and the analytical model where the displacement of IGA model is normalized as:

V=π4kw192QαTL2 (46)

images

Figure 6: Time histories of the non-dimensional deflection at z=L/2 for B=1

It can be found that IGA results are in good agreement with the analytical results even with a very coarse mesh. The good performance of IGA model reveals the superior approximation properties of the high continuous NURBS basis functions which is in contrast to the traditional finite elements.

4.2 A Solar Array of Hubble Space Telescope (HST)

In this example, thermally induced vibrations of the Hubble Space Telescope (HST) solar arrays [8] are investigated. The solar array model consists of two booms, a solar blanket and a spreader bar as shown in Fig. 7. The material properties of the boom model are given in Tab. 1. Besides, the mass of the spreader bar is Ms=1.734 kg and the mass density of the blanket is σsb=1.589 kg/m2.

images

Figure 7: Hubble space telescope solar array model

images

We first studied the thermal behaviors of the boom. The cross section of the boom is a thin-walled annulus with the wall thickness h=2.35×104 m. Due to symmetry, only one half of the annulus is modeled, cf. Fig. 7. The initial temperature of the boom is set to be T0=290 K. For isogeometric thermal analysis, 18 quadratic 1D NURBS elements presented in 3.1 are used, where the exact geometry descriptions of IGA are fully exploited. For comparisons, ABAQUS solutions [52] obtained with 4-node linear heat transfer quadrilateral element DC2D4 are presented, where 1 element and 80 elements are used in the thickness and circumferential directions, respectively. The time histories of the temperature fields of the boom at the locations θ=0 and θ=180 are shown in Fig. 8. In addition, the distributions of the temperature fields around the half circle at different times are shown in Fig. 9. It can be observed from Figs. 8 and 9 that, IGA results agree very well with the ABAQUS reference solutions. The convergence properties of the IGA and ABAQUS models are studied in Tab. 2. It can be found that, IGA results converge faster than ABAQUS results, where only 20 elements are need to reach the converged solutions. While for ABAQUS model, 80 elements are need to reach the converged solutions. This observation confirms the superior properties of isogeometric analysis especially for curved geometries.

images

Figure 8: Time histories of the temperature T at two locations: θ=0, θ=180

images

Figure 9: The distributions of the temperature fields along the circumference of the cross section at different times

images

Based on the above temperature fields, the time histories of the thermal moment can be obtained from Eq. (31) which are shown in Fig. 10. For comparisons, reference solutions taken from Thornton et al. [8] are also shown in Fig. 10. It can be found that, the IGA results predict slightly higher thermal moment compared to the reference solutions (blue curves). This is mainly due to the assumptions adopted in the analytical thermal analysis model given in [8]. The temperature field T(θ,t) in [8] is approximated as the sum of an average temperature T¯(t) and a perturbation temperature Tm(t)cosθ:

T(θ,t)=T¯(t)+Tm(t)cosθ (47)

where the amplitude of the perturbation temperature is assumed to be small in comparison to the average temperature, that is:

Tm/T¯<1 (48)

In addition, the incident heat flux distribution is expanded as a Fourier series neglecting higher order terms as:

δcosθ=(1/π)+12cosθ (49)

Based on the above assumptions, the governing Eq. (9) can be decoupled into two ordinary differential equations:

dT¯dt+σερchT¯4=1παsS0ρch (50)

dTmdt+(kρcR2+4σερchT¯3)Tm=12αsS0ρch (51)

The average temperature in steady-state can be obtained from Eq. (50) at the radiation equilibrium, that is:

T¯ss=(1παsS0σε)14 (52)

When solving for Tm from Eq. (51) a further assumption of T¯=T¯ss is adopted in [8], which simplifies Eq. (51) to a linear differential equation. We think the assumptions adopted in [21] plays an important role for the differences observed in Fig. 10. As mentioned by Thornton [53], the analytical solutions of thermal moment show great differences compared to the finite element solutions especially for lower values of parameter k/(ρcR2). Therefore, a different assumption of T¯=T0 is tested for this problem where the results show slightly difference to the assumption of T¯=T¯ss and agree quite well with the IGA results, cf. Fig. 10. The above observation indicates that the assumption of T¯ is crucial for the analytical method and has to be well defined for specific problems.

images

Figure 10: Time histories of the thermal moment

To proceed the thermally induced vibration analysis, a one-dimensional solar array model is adopted in this example, where the boom and the blanket are coupled at the tip, and where the spreader bar is simplified as a mass point at the tip, cf. Fig. 11. In the HST solar array model, the booms are pre-stressed to maintain the shape of the blanket, therefore, we modify the governing equations of the beam in Eq. (30) by adding an additional term P2w/z2 as:

EI4wz4+2MTz2+ρA2wt2+P2wz2=0 (53)

where P is the boom axial compressive force.

images

Figure 11: Illustration of the coupled 1D HST solar array model

The governing equations of the solar blanket is formulated as:

Fz2wsbz2=σsb2wsbt2 (54)

where wsb is the deflection of the solar blanket, and where Fz is the tension of the solar blanket which can be written as Fz=P/b.

The governing equations of the spreader bar is presented as:

2Vy(L,t)+bbFzwsbz(L,t)dx+Msd2wsdt2=0 (55)

where ws is the deflection of the spread bar, and where Vy is the shear force of boom defined as:

Vy=EI3wbz3PwbzMTz (56)

For isogeometric vibration analysis, both the boom and the solar blanket are discretized with 8 cubic elements, which are the same as numerical Example 1, since the major different between Examples 1 and 2 is thermal analysis. The discretized form of the governing equations of the HST solar array reads:

Msa2wsat2+2ςω0Msawsat+Ksawsa=Fsa (57)

where ς and ω0 are the damping ratio and the first mode natural frequency, respectively, where Msa and Ksa are the mass and stiffness matrices of the coupled system, respectively, where Fsa and wsa are the external force and displacement unknown vectors of the coupled system, respectively.

Various influential factors, e.g., coupling scheme, incident angle and damping ratios, are studied to investigate the thermally induced vibrations of the HST solar array. Fig. 12 shows the tip displacements of the boom for the uncoupled HST model under different pre-stress P, where the critical buckling load Pcr is equal to 48.3 N. In order to evaluate the influences of the analytical thermal analysis on the dynamic behaviors of the solar array, two IGA models are studied, namely the “full-IGA” model and the “semi-IGA” model. In “full-IGA” model, both thermal and structural problems are solved using isogeometric method, while in “semi-IGA” model, only structural problems are solved using isogeometric method, the temperature field of the boom is computed from the analytical method of [8]. It can be found that, “semi-IGA” model agrees very well with the analytical results [8], while “full-IGA” model predicts slightly larger displacements than the analytical results. This discrepancy is mainly due to the approximations introduced in the analytical model for thermal analysis which is proved to be less accurate than IGA models.

images

Figure 12: Thermally induced vibrations of the uncoupled HST solar array under different axial pressures

Figs. 13 and 14 show the time histories of the boom’s tip displacement obtained from the coupled scheme presented in 3.3, where the axial pre-stress is set to be P=14.75 N. It can be found that, thermal flutter of the boom is quite sensitive to the damping ratios and incident angles. For smaller damping ratio and larger incident angle, e.g., ς=0.0001 and θ=80, thermal flutter has been triggered, cf. Fig. 14. While for ς=0.01 and θ=5, thermal flutter can be suppressed, cf. Fig. 13. In addition, the coupling scheme plays an important role in the thermal flutter analysis of the HST solar array, since no such phenomenon can be observed when using the uncoupled scheme, cf. Fig. 12. We note that, when the incident angle is small (θ=5), semi-IGA model agrees very well with the analytical model, while full-IGA model predicts slightly larger displacements which is mainly due to the approximations introduced in the analytical model for thermal analysis. When the incident angle increases to a larger value of θ=80, the discrepancy between full-IGA and analytical model vanishes, which might due to the smaller thermal moments generated in the boom.

images

Figure 13: Stable deflection history for HST solar array from coupled isogeometric thermal-structural analysis

images

Figure 14: Unstable deflection history for HST solar array from coupled isogeometric thermal-structural analysis

4.3 Thin-Walled Lenticular Tubes with Different Cross Sections

The lenticular tubes are widely used in solar sail and antenna structures due to its high efficiency in folding/unfolding operations [54]. In this example, thermally induced vibrations of thin-walled lenticular tubes with different cross sections are studied. The length, thickness and material properties of the tube are chosen the same as the HST solar array problem presented in Section 4.2. The typical cross sectional shapes of a lenticular tube are shown in Fig. 15a, where β is the circular angle which controls the shape of the cross section, where Lf and Ln are the length of the flange and non-flange part, respectively, where h is the thickness of the tube. The thickness of the flange part is 2h and we assume a constant temperature distribution through the thickness of the flange. Diffusive radiation condition is assumed on the outer surface of the tube, while on the internal surface, radiation is neglected. Besides, the initial temperature of the tube is set to be 290 K.

images

Figure 15: Thin-walled lenticular tubes with different cross sections: (a) geometrical parameters of the cross section, (b) coupling constraints between the flange and non-flange parts, where five different shapes of the cross section, namely A (β=20), B (β=30), C (β=45), D (β=60) and E (β=80) are proposed

Five different cross sectional shapes are proposed with the parameter β chosen as β=20, 30, 45, 60, 80, respectively, cf. Fig. 15b. The length of the non-flange part Ln is kept constant so that the solar energy absorbed is the same for the five cross sections. Due to the symmetric boundary conditions, only one half of the tube’s cross section is modeled, and two NURBS patches are used in the thermal analysis of the tube, cf. Fig. 15b. We note that, the complex shapes of the cross section can be modeled exactly using NURBS descriptions which is an intrinsic advantage compared to the traditional finite element method. Coupling constraints need to be applied on the interface of the flange and non-flange part to properly transfer heat fluxes, cf., Fig. 15b. In this paper, we use penalty approach [55] to enforce the coupling constraints T(1)=T(2), where the superscripts ()(1) and ()(2) denote the non-flange and the flange parts, respectively. A moderate penalty value of 104 is chosen while not destroying the conditions of the system matrix. For isogeometric discretization, 19 quadratic NURBS elements are used for thermal analysis while 8 cubic NURBS elements are used for structural analysis.

Fig. 16 shows the temperature distributions along the circular direction of tubes at the time t=2000 s. It can be found that with the increase of β, the average temperature of the tube decreases, however, the difference between the highest and the lowest temperatures increases from 10.18 K for β=20 to 13.88 K for β=80. Fig. 17 shows the time histories of the thermal moments for five different cross sections, it can be found that, thermal moment increases monotonously with the parameter β, which is in contrast to the temperature results in Fig. 16. The main reason behind is, larger values of β increases the height of the cross section which, in consequence, increases the thermal moment of the tube, cf. Eq. (31). We then studied the thermally induced vibrations of the tube with different cross sections in Fig. 18, where larger amplitudes and displacements of the tube are found for smaller values of β. This is simply because, the tube with smaller β has lower moment of inertia which is a dominating factor for the dynamic behaviors of the tube.

images

Figure 16: Temperature distributions along the circular direction of tubes for five cross sections at t=2000 s

images

Figure 17: Time histories of thermal moment with different cross sections

images

Figure 18: Time histories of displacement at the tip for different cross sections

We further proposed two new configurations of the lenticular tube where ribs with a built-in crease are added between the upper and lower part of the cross section, cf. Fig. 19. The ribs are slightly tilted and acting as a hinge to facilitate the folding/unfolding operations. We use two straight NURBS lines to model a rib with a built-in crease. These two lines are C0-continuous at the connection point which is perturbed with a predefined distance, cf. Fig. 19. The parameter β is set to be β=60. For thermal analysis, 19 quadratic NURBS elements are used for the non-flange part and each rib is discretized with four quadratic NURBS elements. The same penalty approach is used to enforce the coupling constraints between the ribs and the non-flange part.

images

Figure 19: Three designs of the lenticular cross section: (I) without ribs, (II) with one rib, (III) with two ribs

The temperature fields of the two new cross section types together with the original design are shown in Fig. 20. ABAQUS solutions obtained with 80 DC2D4 heat transfer elements are presented. It can be found that IGA results agree very well with the ABAQUS reference solutions which demonstrates the accuracy and superior properties of isogeometric analysis. Besides, the temperature gradient of the cross section decreased significantly by adding one and two ribs. Fig. 21 shows the thermal moments of the three designs at the steady state where a reduction of 52.2% and 71.2% of the thermal moment for one and two ribs configurations compared to the original design can be achieved. Additionally, thermally induced vibrations of the newly designed tubes are shown in Fig. 22. Compared to the original design, the maximum displacements of the new designs are reduced by 51.9% and 70.8% for the one and two ribs configurations, respectively. It is worth to mention that, we do not consider the influence of the added ribs on the moment of inertia which results in a conservative prediction of the thermal vibrations.

images

Figure 20: Temperature distributions of the three designs at t = 2000 s

images

Figure 21: Time histories of thermal moment for the three cross sections

images

Figure 22: Time histories of the displacement at the tip for the three cross sections

5  Summary and Conclusions

The flexible appendages of spacecraft are likely to suffer from thermally induced vibrations during day-night or night-day transitions due to the suddenly changed heat load. In this paper, we proposed an isogeometric analysis (IGA) framework for the thermally induced vibration analysis of beam structures where the geometric exactness and higher order continuity of the NURBS basis are fully exploited. For thermal analysis, we developed a one-dimensional curved isogeometric element including radiation heat dissipation, in particular, arbitrary shaped cross sections can be modeled exactly using NURBS geometries. For structural dynamic analysis, the rotation-free Euler–Bernoulli beam elements are used where the C1-continuity requirements can be naturally fulfilled by the NURBS basis functions which is difficult for traditional Lagrange polynomial based finite element method. In addition, we presented a sequential coupling scheme where the thermal and structural responses are solved iteratively.

We studied the accuracy and efficiency of the proposed method with several numerical examples including a benchmark example, a solar array of Hubble Space Telescope (HST) and thin-walled lenticular tubes with different cross sections. Either benchmark or more sophisticated examples show good agreements with the reference solutions. In the HST solar array example, isogeometric thermal element converges much faster than ABAQUS reference solutions due to the higher order, higher continuous and the geometric exactness of the NURBS basis. It is found that, thermal flutters can be well predicted when the thermal-structural coupling is considered, and it is very sensitive to the damping ratios and incident angle. Besides, the slight discrepancy between the IGA and analytical results might due to the approximations of the temperature fields introduced in the analytical methods. In the last example, the influence of cross-sectional shapes on the thermally induced vibrations of the lenticular tube are studied where the modeling of complex cross-sectional shapes is facilitated by the NURBS geometries. The comparisons of the thermal results obtained with IGA and ABAQUS confirm the advantages of the proposed analysis framework. Additionally, two new cross-sectional shapes of the lenticular tube are proposed, numerical results reveal that thermally induced vibrations can be effectively suppressed. It can be concluded that, the proposed framework is suitable for the thermally induced vibration analysis of space structures.

Acknowledgement: The authors are grateful for the important and valuable comments of anonymous reviewers.

Funding Statement: Y. Guo would like to thank the National Natural Science Foundation of China (Grant No. 11972187) and Priority Academic Program Development of Jiangsu Higher Education Institutions for their support.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

  1. Thornton, E. A., Chini, G. P., & Gulik, D. W. (1995). Thermally induced vibrations of a self-shadowed split-blanket solar array. Journal of Spacecraft and Rockets, 32(2), 302-311. [Google Scholar] [CrossRef]
  2. Yu, Y. (1969). Thermally induced vibration and flutter of a flexible boom. Journal of Spacecraft and Rockets, 6(8), 902-910. [Google Scholar] [CrossRef]
  3. Foster, C. L., Tinker, M. L., Nurre, G. S., & Till, W. A. (1995). Solar-array-induced disturbance of the hubble space telescope pointing system. Journal of Spacecraft and Rockets, 32(4), 634-644. [Google Scholar] [CrossRef]
  4. Nurre, G. S., Sharkey, J. P., Waites, H. B. (1991). Initial performance improvements due to design modifications for the pointing control system on the hubble space telescope. Proceedings of the Annual Rocky Mountain Guidance and Control Conference, pp. 493–511, Colorado.
  5. Boley, B. A. (1956). Thermally induced vibrations of beams. Journal of the Aeronautical Sciences, 23(2), 179-181. [Google Scholar] [CrossRef]
  6. Jones, J. P. (1966). Thermoelastic vibrations of a beam. The Journal of the Acoustical Society of America, 39(3), 542-548. [Google Scholar] [CrossRef]
  7. Kraus, H. (1966). Thermally induced vibrations of thin nonshallow spherical shells. AIAA Journal, 4(3), 500-505. [Google Scholar] [CrossRef]
  8. Thornton, E. A., & Kim, Y. A. (1993). Thermally induced bending vibrations of a flexible rolled-up solar array. Journal of Spacecraft and Rockets, 30(4), 438-448. [Google Scholar] [CrossRef]
  9. Beam, R. M. (1969). On the phenomenon of thermoelastic instability (thermal flutter) of booms with open cross section. Washington DC, USA: National Aeronautics and Space Administration.
  10. Rimrott, F. P. J., & Abdel-Sayed, R. (1976). Flexural thermal flutter under laboratory conditions. Transactions of the Canadian Society for Mechanical Engineering, 4(4), 189-196. [Google Scholar] [CrossRef]
  11. Su, X. M., Zhang, J. H., Wang, J., Bi, Y. Q., & Qie, D. F. (2015). Experimental investigation of the thermally induced vibration of a space boom section. Science China Physics, Mechanics and Astronomy, 58(4), 1-9. [Google Scholar] [CrossRef]
  12. Murozono, M., & Thornton, E. A. (1998). Buckling and quasistatic thermal-structural response of asymmetric rolled-up solar array. Journal of Spacecraft and Rockets, 35(2), 147-155. [Google Scholar] [CrossRef]
  13. Mason, J. B. (1968). Analysis of thermally induced structural vibrations by finite element techniques. NASA Technical Report.
  14. Chen, C. S., Chen, C. W., Chen, W. R., & Chang, Y. C. (2013). Thermally induced vibration and stability of laminated composite plates with temperature-dependent properties. Meccanica, 48(9), 2311-2323. [Google Scholar] [CrossRef]
  15. Li, J. L., & Yan, S. Z. (2014). Thermally induced vibration of composite solar array with honeycomb panels in low earth orbit. Applied Thermal Engineering, 71(1), 419-432. [Google Scholar] [CrossRef]
  16. Kumar, R., Mishra, B. K., & Jain, S. C. (2008). Thermally induced vibration control of cylindrical shell using piezoelectric sensor and actuator. The International Journal of Advanced Manufacturing Technology, 38(5–6), 551-562. [Google Scholar] [CrossRef]
  17. Chen, B. S., Gu, Y. X., Zhang, H. W., & Zhao, G. Z. (2003). Structural design optimization on thermally induced vibrations. International Journal for Numerical Methods in Engineering, 58, 1187-1212. [Google Scholar] [CrossRef]
  18. Fan, L. J., & Xiang, Z. H. (2015). Suppressing the thermally induced vibration of large-scale space structures via structural optimization. Journal of Thermal Stresses, 38(1), 1-21. [Google Scholar] [CrossRef]
  19. Xue, M. D., Duan, J., & Xiang, Z. H. (2007). Thermally-induced bending-torsion coupling vibration of large scale space structures. Computational Mechanics, 40(4), 707-723. [Google Scholar] [CrossRef]
  20. Li, W., Xiang, Z. H., Chen, L. J., & Xue, M. D. (2007). Thermal flutter analysis of large-scale space structures based on finite element method. International Journal for Numerical Methods in Engineering, 69, 887-907. [Google Scholar] [CrossRef]
  21. Liu, L., Sun, S. P., Cao, D. Q., & Liu, X. Y. (2019). Thermal-structural analysis for flexible spacecraft with single or double solar panels: A comparison study. Acta Astronautica, 154, 33-43. [Google Scholar] [CrossRef]
  22. Shen, Z. X., Tian, Q., Liu, X. N., & Hu, G. K. (2013). Thermally induced vibrations of flexible beams using absolute nodal coordinate formulation. Aerospace Science and Technology, 29(1), 386-393. [Google Scholar] [CrossRef]
  23. Shen, Z. X., Li, H. J., Liu, X. N., & Hu, G. K. (2017). Thermal shock induced dynamics of a spacecraft with a flexible deploying boom. Acta Astronautica, 141(4), 123-131. [Google Scholar] [CrossRef]
  24. Hughes, T. J. R., Cottrell, J. A., & Bazilevs, Y. (2005). Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194(39–41), 4135-4195. [Google Scholar] [CrossRef]
  25. Bazilevs, Y., Hsu, M. C., Kiendl, W., Wüchner, J. R., & Bletzinger, K. U. (2011). 3D simulation of wind turbine rotors at full scale. Part II-fluid-structure interaction modeling with composite blades. International Journal for Numerical Methods in Fluids, 65, 236-253. [Google Scholar] [CrossRef]
  26. de Lorenzis, L., Temizer, İ., Wriggers, P., & Zavarise, G. (2011). A large deformation frictional contact formulation using NURBS-based isogeometric analysis. International Journal for Numerical Methods in Engineering, 87(13), 1278-1300. [Google Scholar] [CrossRef]
  27. Lezgy-Nazargah, M. (2015). Fully coupled thermo-mechanical analysis of bi-directional FGM beams using NURBS isogeometric finite element approach. Aerospace Science and Technology, 45(6), 154-164. [Google Scholar] [CrossRef]
  28. Zhang, G. D., Alberdi, R., & Khandelwal, K. (2016). Analysis of three-dimensional curved beams using isogeometric approach. Engineering Structures, 117, 560-574. [Google Scholar] [CrossRef]
  29. Guo, Y. (2017). Global-local model coupling for composite shell structures in the framework of isogeometric analysis. Composite Structures, 176(14), 617-629. [Google Scholar] [CrossRef]
  30. Guo, Y., Do, H., & Ruess, M. (2019). Isogeometric stability analysis of thin shells: from simple geometries to engineering models. International Journal for Numerical Methods in Engineering, 118(8), 433-458. [Google Scholar] [CrossRef]
  31. van Do, V. N., & Lee, C. H. (2019). Quasi-3D isogeometric buckling analysis method for advanced composite plates in thermal environments. Aerospace Science and Technology, 92(83), 34-54. [Google Scholar] [CrossRef]
  32. Nguyen-Thanh, N., Li, W., & Zhou, K. (2018). Static and free-vibration analyses of cracks in thin-shell structures based on an isogeometric-meshfree coupling approach. Computational Mechanics, 62(6), 1287-1309. [Google Scholar] [CrossRef]
  33. Nagy, A. P., Abdalla, M. M., & Gürdal, Z. (2010). Isogeometric sizing and shape optimisation of beam structures. Computer Methods in Applied Mechanics and Engineering, 199(17–20), 1216-1230. [Google Scholar] [CrossRef]
  34. Weeger, O., Wever, U., & Simeon, B. (2013). Isogeometric analysis of nonlinear Euler–Bernoulli beam vibrations. Nonlinear Dynamics, 72(4), 813-835. [Google Scholar] [CrossRef]
  35. Bauer, A. M., Breitenberger, M., Philipp, P., Wüchner, R., & Bletzinger, K. U. (2016). Nonlinear isogeometric spatial Bernoulli beam. Computer Methods in Applied Mechanics and Engineering, 303, 101-127. [Google Scholar] [CrossRef]
  36. Kiendl, J., Bletzinger, K. U., Linhard, J., & Wüchner, R. (2009). Isogeometric shell analysis with Kirchhoff–Love elements. Computer Methods in Applied Mechanics and Engineering, 198(49–52), 3902-3914. [Google Scholar] [CrossRef]
  37. Oesterle, B., Ramm, E., & Bischoff, M. (2016). A shear deformable, rotation-free isogeometric shell formulation. Computer Methods in Applied Mechanics and Engineering, 307, 235-255. [Google Scholar] [CrossRef]
  38. Duong, T. X., Roohbakhshan, F., & Sauer, R. A. (2017). A new rotation-free isogeometric thin shell formulation and a corresponding continuity constraint for patch boundaries. Computer Methods in Applied Mechanics and Engineering, 316(1–3), 43-83. [Google Scholar] [CrossRef]
  39. Guo, Y., Heller, J., Hughes, T. J. R., Ruess, M., & Schillinger, D. (2018). Variationally consistent isogeometric analysis of trimmed thin shells at finite deformations, based on the STEP exchange format. Computer Methods in Applied Mechanics and Engineering, 336(49–52), 39-79. [Google Scholar] [CrossRef]
  40. Guo, Y., Zou, Z., & Ruess, M. (2021). Isogeometric multi-patch analyses for mixed thin shells in the framework of non-linear elasticity. Computer Methods in Applied Mechanics and Engineering, 380(49–52), 113771. [Google Scholar] [CrossRef]
  41. Li, W., Nguyen-Thanh, N., & Zhou, K. (2018). Geometrically nonlinear analysis of thin-shell structures based on an isogeometric-meshfree coupling approach. Computer Methods in Applied Mechanics and Engineering, 336(9), 111-134. [Google Scholar] [CrossRef]
  42. Duvigneau, R. (2009). An introduction to isogeometric analysis with application to thermal conduction. INRIA Technical Report.
  43. Chen, T., Mo, R., & Zhang, X. (2011). Isogeometric analysis for transient heat conduction of solid medium. Computer Integrated Manufacturing Systems, 17(9), 1988-1996. [Google Scholar] [CrossRef]
  44. Cottrell, J. A., Reali, A., Bazilevs, Y., & Hughes, T. J. R. (2006). Isogeometric analysis of structural vibrations. Computer Methods in Applied Mechanics and Engineering, 195(41–43), 5257-5296. [Google Scholar] [CrossRef]
  45. Thai, C. H., Nguyen-Xuan, H., Nguyen-Thanh, N., Le, T. H., & Nguyen-Thoi, T. (2012). Static, free vibration, and buckling analysis of laminated composite Reissner–Mindlin plates using NURBS-based isogeometric approach. International Journal for Numerical Methods in Engineering, 91(6), 571-603. [Google Scholar] [CrossRef]
  46. Zhang, X. K., Xia, Y., Hu, Q. Y., & Hu, P. (2017). Efficient isogeometric formulation for vibration analysis of complex spatial beam structures. European Journal of Mechanics-A/Solids, 66(SI), 212-231. [Google Scholar] [CrossRef]
  47. Maurin, F., Dedé, L., & Spadoni, A. (2015). Isogeometric rotation-free analysis of planar extensible-elastica for static and dynamic applications. Nonlinear Dynamics, 81(1–2), 77-96. [Google Scholar] [CrossRef]
  48. Sickinger, C., Herbeck, L., & Breitbach, E. (2006). Structural engineering on deployable CFRP booms for a solar propelled sailcraft. Acta Astronautica, 58(4), 185-196. [Google Scholar] [CrossRef]
  49. Piegl, L., Tiller, W. (1996). The NURBS book. Germany: Springer Science & Business Media.
  50. Cottrell, J. A., Hughes, T. J. R., Bazilevs, Y. (2009). Isogeometric analysis: Towards integration of CAD and FEM. USA: John Wiley & Sons.
  51. Huebner, K. H., Dewhirst, D. L., Smith, D. E., Byrom, T. G. (2001). The finite element method for engineers. USA: John Wiley & Sons.
  52. Hibbit, D., Karlsson, B., Sorensen, P. (2007). ABAQUS/Standard analysis user’s manual. USA: Hibbit, Karlsson, Sorensen Inc.
  53. Thornton, E. A. (1996). Thermal structures for aerospace applications. USA: American Institute of Aeronautics and Astronautics.
  54. Block, J., Straubel, M., & Wiedemann, M. (2011). Ultralight deployable booms for solar sails and other large gossamer structures in space. Acta Astronautica, 68(7–8), 984-992. [Google Scholar] [CrossRef]
  55. Guo, Y., & Tzou, H. (2018). Ultraviolet-activated frequency control of beams and plates based on isogeometric analysis. Journal of Vibration and Acoustics, 140(3), 879. [Google Scholar] [CrossRef]
images This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.