Dynamic Analysis of Spatial Truss Structures Including Sliding Joint Based on the Geometrically Exact Beam Theory and Isogeometric Analysis

With increasing of the size of spatial truss structures, the beam component will be subjected to the overall motion with large deformation. Based on the local frame approach and the geometrically exact beam theory, a beam finite element, which can effectively reduce the rotational nonlinearity and is appropriate for finite motion and deformation issues, is developed. Dynamic equations are derived in the Lie group framework. To obtain the symmetric Jacobian matrix of internal forces, the linearization operation is conducted based on the previously converged configuration. The iteration matrix corresponding to the rotational parameters, including the Jacobian matrix of inertial and internal forces in the initial configuration, can be maintained in the simulation, which drastically improves the computational efficiency. Based on the Lagrangian multiplier method, the constraint equation and its Jacobian matrix of sliding joint are derived. Furthermore, the isogeometric analysis (IGA) based on the non-uniform rational B-splines (NURBS) basis functions, is adopted to interpolate the displacement and rotation fields separately. Finally, three dynamic numerical examples including a deployment dynamic analysis of spatial truss structure are conducted to verify the availability and the applicability of the proposed formulation.


Introduction
The deployable spatial truss structure is a typical flexible multibody system consisting of beam members. With the increase of the size, the deployable spatial truss structure usually undergoes finite deformation when refer to finite displacement and rotation. The most common modeling methods in flexible multibody system for beams are the Absolute Nodal Coordinate Formulation (ANCF) and the Geometrically Exact Beam Theory (GEBT). Both approaches can provide the solution of finite deformation and finite motion. The ANCF, proposed by Shabana [1,2], is one of the methods aimed at multibody system dynamics areas recent years. The method uses slopes instead of infinitesimal rotations to simplify the parameter updating algorithm, but at the cost of introducing more degrees of by the generalized-α method [30,31]. The sliding joint is modeled using the Lagrangian multiplier method. The constraints equations relative to the sliding joint and its Jacobian matrix are also given.
The paper is organized as follows: The basic theory of the geometrically exact beam theory including beam kinematics, static equations, dynamic equations, and their linearization are introduced in Section 2. Finite element strategies are employed in Section 3. In Section 4, the modeling of sliding joint is conducted. In Section 5, two dynamic numerical examples and a deployment dynamic analysis of spatial truss structure are implemented to verify the availability and the applicability of the formulation proposed in this paper. Conclusions are presented in Section 6.

Beam Kinematics
In this section, a straight beam element with a square cross section was used as an example, and this method can be easily extended to curved beam elements with arbitrary cross sections. A reference configuration B 0 and a current configuration B are displayed in Figure 1. For an arbitrary cross-section, the reference point of the local frame is assumed to coincide with the centroid O(o 2 , o 3 ). A unit orthogonal basis system A = {A 1 , A 2 , A 3 } along with the material coordinates {ξ 1 , ξ 2 , ξ 3 } is introduced. The axis of the beam is initially along A 1 , with the length of beam l and the arc-length parameter S = ξ 1 ∈ [0, l] of the spatial curve. The orientation of the cross-section is described by the basis vectors {A 2 , A 3 }. The local frame a = {a 1 , a 2 , a 3 } is attached to the cross-section in the current configuration, which is characterized by the time parameter t. For convenience, it is assumed that the material frame A i and inertial frame e i initially coincide, and thus the relationship between the basis vectors A i and a i is as follows: where the operator R is the cross-section orthogonal rotation tensor with respect to inertial frame. In general, the rotation tensor is defined from the exponential map and can be expressed directly by Rodrigues formula [32]: R = exp(Θ) = I 3×3 + a 1 (Θ)Θ + a 2 (Θ)ΘΘ (2) with Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 20 The paper is organized as follows: The basic theory of the geometrically exact beam theory including beam kinematics, static equations, dynamic equations, and their linearization are introduced in Section 2. Finite element strategies are employed in Section 3. In Section 4, the modeling of sliding joint is conducted. In Section 5, two dynamic numerical examples and a deployment dynamic analysis of spatial truss structure are implemented to verify the availability and the applicability of the formulation proposed in this paper. Conclusions are presented in Section 6.

Beam Kinematics
In this section, a straight beam element with a square cross section was used as an example, and this method can be easily extended to curved beam elements with arbitrary cross sections. A reference configuration 0 and a current configuration are displayed in Figure 1. For an arbitrary cross-section, the reference point of the local frame is assumed to coincide with the centroid along with the material coordinates { , , }  a a a a is attached to the cross-section in the current configuration, which is characterized by the time parameter t. For convenience, it is assumed that the material frame i A and inertial frame i e initially coincide, and thus the relationship between the basis vectors and i a is as follows: where the operator R is the cross-section orthogonal rotation tensor with respect to inertial frame. In general, the rotation tensor is defined from the exponential map and can be expressed directly by Rodrigues formula [32]: the skew-symmetric matrix of the rotational vector Θ and Θ means the magnitude of the rotational vector. On the contrary, the inverse map of the exponential map is the logarithmic map [33], giving: with θ = arccos 1 2 (trace(R) − 1) and the absolute value |θ| < π. An arbitrary point p on the cross-section in the current configuration can be expressed by: in which y p = ξ 2 A 2 + ξ 3 A 3 , and is defined from the reference configuration since the cross-section is assumed to remain undeformed.

Beam Statics and Its Linearization
According to continuum mechanics, the Green-Lagrangian strain tensor E(ξ 1 , ξ 2 , ξ 3 ) is defined as: in which F k j is the deformation gradient, and δ ij is the Kronecker delta. The derivatives of the displacement field in Equation (5) with respect to ξ i (i = 1, 2, 3) are given by: where (·) denotes the derivative relative to arc-length parameter S.
In the classical geometrically exact beam theory, the strain measures ε can be represented as: For an initially straight beam, the strain measures ε can be thus written as: in which γ = [γ 1 , γ 2 , γ 3 ] T quantifies the beam axis extension and two shearing; κ = [κ 1 , κ 2 , κ 3 ] T measures its twist and two curvatures in local frame. The tangent operator can be obtained from [10], giving: with the coefficients: Appl. Sci. 2020, 10, 1231

of 19
Consequently, based on the rigid perimeter and small strain assumptions, the non-vanishing terms in Equation (6) can be given as: Assuming a constitutive law of isotropic linear elastic materials, the second Piola-Kirchhoff stress tensor S = [S 11 , S 12 , S 13 ] T refers to: in which E, and G are the elasticity modulus and shear modulus, respectively. Based on the Green-Lagrange strain tensor, the internal virtual work can be computed as: withĒ = ε and the elastic matrix: where τ is the shear coefficient. The notation A in Equation (15) denotes the area of the cross-section and four constant section properties are calculated from the following formulas: The stress measures can be given by: in which the stress resultant N = [N 1 , N 2 , N 3 ] T represents one normal force and two shear forces in local frame; the stress couple M = [M 1 , M 2 , M 3 ] T includes one torsional moment and two bending moments in local frame. According to Equation (9) the variation of the strain measures can be written as: with δh = [δx T , δΘ T ] T and the strain operator: Appl. Sci. 2020, 10, 1231 6 of 19 In each iterative step, the computed rotation increments lie in the tangent space of the current configuration. To obtain a symmetric tangent stiffness, as suggested by [10], the rotation variation in local frame is transformed in terms of the rotation increment ϕ with respect to previously converged configuration with the relationship: δΘ = Tδϕ (20) Then the Equation (18) can be further expressed as: with the generalized displacement variations δh * = [δx T , δϕ T ] T and the strain operator: Then the internal virtual work can be further written as: Thus, the tangent stiffness matrix can be derived by linearization of Equation (23), giving: (24) in which the first term is the material part and the second term is the geometric part. By introducing an operator Ξ, the continuous geometric stiffness matrix can be further obtained, giving: where H is computed by the formula: with the coefficients, For arbitrary 3 × 1 vector u, the operator Ξ can be expressed as: The matrices Π in the material part are symmetric. The Υ in the geometric part obtained here shows non-obvious symmetry in its formula structure, although it is numerically symmetric. In the case of small deformation, the geometric stiffness is unnecessary in general. However, numerical examples show that the application of geometric stiffness can improve the convergence of Newton iterations for cases of both small deformation and finite deformation.

Beam Dynamics and Its Linearization
The weak form of dynamic equilibrium equation can be expressed as: with the internal virtual work δW int , external virtual work δW ext , and inertial virtual work δW ine . Denoting ρ as the mass density of the material, the inertial virtual work δW ine is given by: where A ρ = Aρ is the density of the cross-section. α and Ω are the angular acceleration and angular velocity of beam cross-section in local frame, respectively. Besides, the notation J denotes the cross-sectional moment of inertia with relative to the centroid, giving: Then the generalized mass matrix can be derived by linearization of the inertial virtual work δW ine , giving: Since the virtual rotational vector δϕ is defined in the local frame, it is free from the influence of rigid rotation, which greatly reduces the rotational nonlinearity. Consequently, the parts of the stiffness matrix and the generalized mass matrix, which corresponds to the rotation parameters, are invariable under arbitrary rigid motion. Thereby, the number of updates required during the computation process decreases sharply.
Similarly, the external virtual work can be computed by: where F andF represent the distributed force and the concentrated force on the centerline of the beam in the inertial frame, respectively. P andP represent the distributed moment and the concentrated moment on the cross-section of the beam in the local frame, respectively.

Spatial Discretization: NURBS Basis Functions
In general, the dynamics equilibrium equations of the geometrically exact beam can be transformed into a set of nonlinear algebraic equations by finite element spatial discretization and time discretization. In the spatial domain, the linearized governing equations described above are discretized by adopting the IGA; namely, both the geometry and the unknown variables are approximated by NURBS [24].
Denote Ψ as a non-decreasing knot vector, giving: with the degree of basis functions ς and the number of control point n. According to the Cox-de Boor recursion formula, a B-spline basis function can be obtained from a linear combination of two B-spline basis functions of degree ς-1, giving: with when ς = 0. Then, a NURBS basis function can be derived through introducing the weight factors w i , giving: in which χ i,ς is the B-spline basis function. Obviously, the NURBS basis function degenerates into the B-spline basis function when all weight factors are equal to one. The position coordinate vector of any point can be calculated from interpolation and written as: where P i is the control point. The continuity of the NURBS basis function is determined by the multiplicity of a knot. Given multiplicity k of a knot, the continuity of the basis function is C ∞ between knots, while it is C ς−k at the knot. Thus, NURBS can accurately describe various types of curves and surfaces, such as circles, hyperbolas, and elliptic surfaces.

Discretization of the Linearized Equations
It is worth noting that separate interpolations of the displacement and rotation fields are conducted to make the geometric representations uniform in the finite element analysis and CAD. Since the displacement field takes values in linear space R 3 while the rotation field in nonlinear manifold SO (3), the interpolations of the displacement and rotation field are conducted in different methods. To satisfy the objectivity of strain measures, the Crisfield and Jelenić interpolation algorithm [12] is adopted, in which the relative rotation is interpolated. The interpolation formulas are given by: where N I is the NURBS basis function corresponding to control point I, and a summation is extended to all the control points of the element. The notation Θ r I is the relative rotation, and R r is the rotation matrix of the reference point, which is typically chosen as the first control point on an element. The strain measures can thus be interpolated as: Hence, the generalized displacement variation and its derivative can be discretized into: where δh * I is the generalized displacement variation of the control point. The strain-displacements matrix B * I is the discrete counterpart of the strains operator B * , as shown in Equation (22), giving: Consequently, the strain measures can be written as: The generalized internal force and inertial force can be discretized from Equations (23) and (30) as: In this work, as a separate interpolation is employed, shear locking exists in 1-order elements and is circumvented by the selective reduced integration. The one-point quadrature strategy is employed to calculate the shear strain, while the full quadrature, i.e., two-point quadrature for 1-order element is used to calculate the membrane and bending strains. Besides, high-order elements will not suffer from the shear locking problem and the full quadrature is employed to obtain all the strains in the Equations (46) and (47).

Time Discretization: Lie Group Generalized-α Method
The dynamics equilibrium equations can be expressed by a set of Differential Algebraic Equations (DAEs), giving: where F ext is the generalized external force; the notation Φ and Φ q denote the constraint equations and its Jacobian matrix relative to the generalized coordinates; and λ is the Lagrange multiplier vector. For time discretization, the generalized-α method is employed to solve the DAEs. The displacement field is handled with the conventional generalized-α method [30] while the Lie group time integrators [31] are adopted to solve the rotation field. The solvers of the rotation field are based on the equations: in which Θ n represents the rotational vector from time t n to time t n+1 ; h = t n+1 − t n means the time step; α m , α f , γ, and β are algorithmic parameters for the generalized-α method. According to the simultaneous Equations (48) and (49), the rotation matrix R n+1 , angular velocity Ω n+1 , and angular acceleration α n+1 at time t n+1 can be obtained through newton iterations.

Sliding Joint
In this work, the sliding joint is implemented using the Lagrange multiplier method. The constraint of sliding joint is composed of three position constraint equations and one constraint equation, which is expressed according to three Lagrange multipliers and one non-generalized degree of freedom that represents the arc length parameter ζ. The complete constraint equations can be written as: where the vectors x A and x B denote the position of joint point on slider and carrier beam, respectively. The constrain Φ in Equation (50) actually leads to a spherical joint, and the Lagrange multiplier λ = [λ 1 , λ 2 , λ 3 ] T is employed to connect with the constraint forces applied to both bodies. When the vector λ is perpendicular to the sliding direction, the constraint force in this direction is released so that the joint point can move along the midline of carrier beam under the frictionless condition.
In the simulation, denotingq T = [q T , ζ] as the generalized coordinates, the dynamics equilibrium equations in Equation (48) would be further written as: The Jacobian matrix of the equilibrium equation can be evaluated as: in which the matrices K int , K ine , and K ext are interpreted as the Jacobian matrix of generalized internal force, generalized inertia force, and generalized external force. The notation (•) ,q denotes the derivative relative to generalized coordinates q. Note that the position vector x B is calculated from interpolation, giving: The subscript I is determined by the interpolation function. It is remarked that for the case of sliding along very flexible beam, the IGA method takes more advantage on accuracy and efficiency when using the C k -continuous NURBS basic function (k 1).

Numerical Examples
Three typical dynamic examples are adopted in this paper. The first one is a high-speed flexible slider-crank with moderate deformation, which tests the proposed method's character of reducing the rotational nonlinearity. The second one is a typical testing case of the sliding joint, which demonstrates that the present formulation is applicable to the analysis of the structures including the sliding joint, and undergoing finite deformation and motion. Finally, the dynamics of a deployable truss structure unit is investigated based on the proposed formulation. The IGA based on the high-order NURBS basis function is applied in all examples.

High-Speed Flexible Slider-Crank
This example is about a flexible slider-crank suffered from a constant moment M, which aims at studying the influence of proposed method on the rotational nonlinearity. As shown in Figure 2, Rod I is connected to the earth and Rod II through a spherical joint at the same time; the rigid slider was also hinged with Rod II by a spherical joint and is restricted to moving along the X-axis. The mass of the slider was 0.1 kg and the moment was M = 500 Nm. Initially, Rod I and Rod II were both horizontally stationary in the X-Y plane, that is, α = 0. Other geometrical and material parameters are shown in Table 1. In the simulation, both rods are discretized into six elements using the 3-order NURBS basic function; the total simulation time was 0.2 s with two different time steps h = 0.0001 s and h = 0.0005 s. The elastic deformation of rod was measured by the distance between the midpoint of rod and the midpoint of the line connecting the two ends of rod. The deformation rate was the ratio between the elastic deformation of the rod and the original length of the rod. Obviously, the deformation of Rod II during the simulation increased gradually with the increase of rotation speed, as displayed in Figure 3. The maximum deformation of Rod II and the maximum rotational speed of Rod I exceeded 5% of the rod length and 320 rad/s, respectively.

High-Speed Flexible Slider-Crank
This example is about a flexible slider-crank suffered from a constant moment M, which aims at studying the influence of proposed method on the rotational nonlinearity. As shown in Figure 2, Rod I is connected to the earth and Rod II through a spherical joint at the same time; the rigid slider was also hinged with Rod II by a spherical joint and is restricted to moving along the X-axis. The mass of the slider was 0.1 kg and the moment was M = 500 Nm. Initially, Rod I and Rod II were both horizontally stationary in the X-Y plane, that is, α = 0. Other geometrical and material parameters are shown in Table 1. In the simulation, both rods are discretized into six elements using the 3-order NURBS basic function; the total simulation time was 0.2 s with two different time steps h = 0.0001 s and h = 0.0005 s. The elastic deformation of rod was measured by the distance between the midpoint of rod and the midpoint of the line connecting the two ends of rod. The deformation rate was the ratio between the elastic deformation of the rod and the original length of the rod. Obviously, the deformation of Rod II during the simulation increased gradually with the increase of rotation speed, as displayed in Figure 3. The maximum deformation of Rod II and the maximum rotational speed of Rod I exceeded 5% of the rod length and 320 rad/s, respectively.   In general, if the Jacobian matrix corresponding to the external force is ignored, then the system iteration matrix consists of a generalized mass matrix, a tangent stiffness matrix, and a  In general, if the Jacobian matrix corresponding to the external force is ignored, then the system iteration matrix consists of a generalized mass matrix, a tangent stiffness matrix, and a Jacobian matrix corresponding to the constraint equation. It is assumed that J θ represents the iterative matrix composed of the generalized mass matrix and the stiffness matrix corresponding to rotation parameters, and the convergence criterion of Newton iteration is set to the incremental value of the generalized coordinate, . The number of times that the iteration matrix J θ is updated (update times) and the total iterations in the simulation for different time discretizations and update conditions of the truss structure unit are given in Table 2. The "Update" means J θ is updated all the time, and the "Frozen" denotes J θ in the reference configuration is maintained for the whole simulation. Obviously, if J θ is updated in each Newton iteration step, the whole simulation needs a total of 1435 Newton iteration steps, and J θ needs to be updated 1434 times. While the iteration matrix J θ is kept for the whole simulation, the total iterations merely needs 1473 times. It is demonstrated that the rigid rotational nonlinearity is reduced dramatically based on the local frame approach, such that the iteration matrix J θ can remain unchanged all the time, especially for the case of moderate deformation rate in this example. Moreover, the rigid rotational nonlinearity can be also reduced with a smaller time step, since the total iterations remains unchanged in the condition of "Frozen", as shown in Table 2. From a computational point of view, however, the smaller the time step, the more iteration steps are needed, which increases the computation cost.

Parameters
First  In general, if the Jacobian matrix corresponding to the external force is ignored, then the system iteration matrix consists of a generalized mass matrix, a tangent stiffness matrix, and a X Y M  Figure 3. The deformation rate of Rod II. Table 2. The number of times that the iteration matrix is updated (update times) and the total iterations in the simulation for different time discretizations and update conditions of the truss structure unit. Update: the iteration matrix is updated all the time. Frozen: the iteration matrix in the reference configuration is maintained for the whole simulation.

Sliding Three-Dimensional Beam with Eccentricity
Two initial straight flexible beams connected through a sliding joint are involved in this example, as shown in Figure 4. Two spherical joints locate at the two ends of the carrier beam AB. A concentrated mass M is fixed at the free end of sliding beam AM. All the geometrical and material parameters are displayed in Table 3. This is a typical test case of sliding joint, which has been proposed by Sugiyama et al. [19] and applied by other researchers [20,23]. The ANCF and GEBT are employed in [19,20] and [23], respectively. In reference [19,23], the carrier beam and sliding beam are both discretized into one element thereby obtaining an un-converged solution. The convergence of solution was investigated in [23].
terminated once the sliding beam reaches point B, which occurred at t = 1.306 s. The deformed configurations during the motion are shown in Figure 5 with the time interval Δt = 0.1 s, while the X-Z and Y-Z trajectories of the mass are plotted in Figure 6. It is shown that the flexible beams have visible deformation during the finite translational and rotational process. The results from the present method have a good agreement with the reference solution from ANCF. It is demonstrated that the present formulation is applicable to the analysis of the structures with sliding joint, which undergo finite deformation and motion.     In the example, the carrier beam and sliding beam are discretized into 16 and 2 elements using the 2-order NURBS basic function. For comparison, the reference solutions from ANCF with 64 and 4 elements were calculated. The calculation was conducted with the time step h = 0.001 s and terminated once the sliding beam reaches point B, which occurred at t = 1.306 s. The deformed configurations during the motion are shown in Figure 5 with the time interval ∆t = 0.1 s, while the X-Z and Y-Z trajectories of the mass are plotted in Figure 6. It is shown that the flexible beams have visible deformation during the finite translational and rotational process. The results from the present method have a good agreement with the reference solution from ANCF. It is demonstrated that the present formulation is applicable to the analysis of the structures with sliding joint, which undergo finite deformation and motion.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 20 In the example, the carrier beam and sliding beam are discretized into 16 and 2 elements using the 2-order NURBS basic function. For comparison, the reference solutions from ANCF with 64 and 4 elements were calculated. The calculation was conducted with the time step h = 0.001 s and terminated once the sliding beam reaches point B, which occurred at t = 1.306 s. The deformed configurations during the motion are shown in Figure 5 with the time interval Δt = 0.1 s, while the X-Z and Y-Z trajectories of the mass are plotted in Figure 6. It is shown that the flexible beams have visible deformation during the finite translational and rotational process. The results from the present method have a good agreement with the reference solution from ANCF. It is demonstrated that the present formulation is applicable to the analysis of the structures with sliding joint, which undergo finite deformation and motion.

Deployable Spatial Truss Structure Unit
The dynamics of a deployable truss structure unit, which can be employed to construct large space systems [34], was investigated in this section. As shown in Figure 7, the structure unit is an assembly composed of three square frames, eight cross-bars, and eight diagonal rods joined to one another. The cross-bars and the diagonal rods on both sides are jointed with square frames with only the rotations around the ξ 2 -axis released, while the above and below diagonal rods are clamped at corresponding points on square frames. During the deployment process, the diagonal rods on both sides can slide between two ends of the carrier rods through a sliding joint. The geometrical and material parameters can be found in Table 4.
the rotations around the 2  -axis released, while the above and below diagonal rods are clamped at corresponding points on square frames. During the deployment process, the diagonal rods on both sides can slide between two ends of the carrier rods through a sliding joint. The geometrical and material parameters can be found in Table 4. In the first test case, the first (left) square frame was fixed initially, and the structure unit was kept with i  = 5° (i = 1, 2, 3, 4). Then, the structure unit was driven with the gravity 0.005 g along X-axis. The calculation was conducted with the time step h = 0.1 s and a single 2-order NURBS finite element per beam. When i  = 90° (i = 1, 2, 3, 4), as shown in Figure 7, the sliding joints were locked in the top end of carrier beams, and the calculation was terminated through a constraint of ζ. For comparison, the reproductive simulation based on MSC.ADAMS and ANSYS were conducted. As suggested by [35,36], the modal neutral files (.mnf files) was firstly obtained from ANSYS and was then imported into ADAMS/FLEX module, which allows us to model flexible parts based on its modal frequency. The beam members were modeled with the specified parameters in Table 4 using Solid45 element in ANSYS. The X-position of point D during the motion calculated from proposed formulation were of good agreement with MSC.ADAMS, as displayed in Figure 8. It was demonstrated that the present formulation was appropriate for the deployment dynamic analysis of the spatial truss structure.   In the first test case, the first (left) square frame was fixed initially, and the structure unit was kept with β i = 5 • (i = 1, 2, 3, 4). Then, the structure unit was driven with the gravity 0.005 g along X-axis. function is the best control strategies in four functions, since it shows the best deployment stability and minimum residual vibration after the deployment.

Conclusions
In order to investigate the dynamics of the deployable spatial truss structure, a beam finite element that can effectively reduce the rotational nonlinearity thereby applicable for finite motion and deformation issues, was proposed based on the geometrically exact beam theory and the local frame approach. The analytical expressions of internal and inertial forces and their Jacobian matrices, including the symmetric stiffness matrix, were strictly derived in this paper. Based on the Lagrangian multiplier method, the constraint equation and its Jacobian matrix of sliding joint were also derived. Using NURBS basic function, the displacement and rotation fields were interpolated separately to be compatible with IGA. The numerical example illustrated that the iteration matrix corresponding to the rotational parameters in the initial configuration, including the Jacobian matrix of inertial and internal forces, could be maintained in the whole simulation, which drastically improved the computational efficiency. Moreover, the sliding joint was also successfully tested. The application example demonstrated that the present method was applicable to the deployable truss structure with sliding joint. The simulation results show that the effect of the cosine function on angle planning was better than the linear function, polynomial function, and cycloid function.
Author Contributions: J.R. and C.L. conceived and designed the project; Z.W., C.L., Z.L., P.X., W.L., and W.S. designed the numerical experiments; Z.W. performed the computations and analyzed the results; Z.L. and P.X. contributed analysis tools; Z.W., J.R. and C.L. wrote the paper.

Conclusions
In order to investigate the dynamics of the deployable spatial truss structure, a beam finite element that can effectively reduce the rotational nonlinearity thereby applicable for finite motion and deformation issues, was proposed based on the geometrically exact beam theory and the local frame approach. The analytical expressions of internal and inertial forces and their Jacobian matrices, including the symmetric stiffness matrix, were strictly derived in this paper. Based on the Lagrangian multiplier method, the constraint equation and its Jacobian matrix of sliding joint were also derived. Using NURBS basic function, the displacement and rotation fields were interpolated separately to be compatible with IGA. The numerical example illustrated that the iteration matrix corresponding to the rotational parameters in the initial configuration, including the Jacobian matrix of inertial and internal forces, could be maintained in the whole simulation, which drastically improved the computational efficiency. Moreover, the sliding joint was also successfully tested. The application example demonstrated that the present method was applicable to the deployable truss structure with sliding joint. The simulation results show that the effect of the cosine function on angle planning was better than the linear function, polynomial function, and cycloid function.
Author Contributions: J.R. and C.L. conceived and designed the project; Z.W., C.L., Z.L., P.X., W.L., and W.S. designed the numerical experiments; Z.W. performed the computations and analyzed the results; Z.L. and P.X. contributed analysis tools; Z.W., J.R. and C.L. wrote the paper. All authors have read and agreed to the published version of the manuscript.