Investigation of strut collision in tensegrity statics and dynamics

International Journal of Solids and Structures 167 (2019) 202–219

Contents lists available at ScienceDirect

International Journal of Solids and Structures journal homepage: www.elsevier.com/locate/ijsolstr

Investigation of strut collision in tensegrity statics and dynamics Ziyun Kan a, Haijun Peng a,∗, Biaoshong Chen a, Xiaohui Xie b, Lining Sun b a

Department of Engineering Mechanics, State Key Laboratory of Structural Analysis for Industrial, Equipment, Dalian University of Technology, Dalian 116024, China b School of Mechanical and Electric Engineering, Soochow University, Suzhou 215021, China

a r t i c l e

i n f o

Article history: Received 24 May 2018 Revised 23 January 2019 Available online 12 March 2019 Keywords: Tensegrity Strut collision Dynamic deployment Multibody system Sliding cables

a b s t r a c t Strut collisions can occur in tensegrity structures and seriously change the general mechanical behaviors of the structures. For many years, there have been few investigations into this topic. This paper systematically addresses the issue of strut collision modeling in tensegrity mechanical analyses. A general and efficient rigid body model for strut collision analysis is presented. A two-strut collision model is formulated using contact force methods. The normal contact force and the tangential friction force are evaluated by using the Lankarani–Nikravesh model and the modified Coulomb’s friction law, respectively. Exact formulations of the generalized force vectors and the related tangent stiffness and tangent damping matrices are analytically derived. Fundamental issues related to collision detection and stiffness coefficient determination, are also addressed. Two representative examples involving both classical and clustered tensegrity structures are presented to demonstrate the effectiveness of the general framework, as well as to make a detailed study into the influence of collisions effects. Outcomes presented in this work will broaden the research discipline on tensegrity systems. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction In recent years, new trends of engineering applications have excited researchers’ interests in new structural style designs. Among many new structural styles, the tensegrity concept may be one of the most promising for new structural designs (Nagase and Skelton, 2015; Sultan, 2009). Although the precise definition of tensegrity structures is subject to debate (Hanaor, 2012), tensegrity structures considered here are composed of struts and cables, and their integrity is maintained by balancing tensional forces in cables and compressive forces in struts (Snelson, 1965). Due to the existence of a large number of lightweight cables, tensegrity structures usually have a high strength-to-mass ratio, which leads to the development of a strong and lightweight structural design (Skelton and Oliveira, 2009). Tensegrity as a concept is particularly attractive for a wide variety of technologically important applications in civil architectures, mechanical devices, robotics and biomechanics (Skelton and Oliveira, 20 09; Sultan, 20 09). Mechanical analyses of tensegrity structures have been conducted for many years (Calladine, 1978; Kebiche et al., 1999; Pellegrino, 2001). The influence of various kinds of loads involved in tensegrity systems has been investigated, such as prestress forces (Murakami, 2001; Tran and Lee, 2010), distributed ∗

Corresponding author. E-mail address: [email protected] (H. Peng).

https://doi.org/10.1016/j.ijsolstr.2019.03.012 0020-7683/© 2019 Elsevier Ltd. All rights reserved.

and concentrated external forces (including gravitational forces) (Faroughi et al., 2015), and recently, thermal loads (Ashwear and Eriksson, 2015). Such studies may be conducted from a static or dynamic point of view and undoubtedly offer a deeper understanding of such structures. Despite the large body of literature addressing mechanical analysis, one type of load that has seldom been examined: the force induced by strut collision. As a result of such a force, the state of structures can change quite rapidly, resulting in discontinuities of system velocities, accompanied by vibrations and load propagations in structural components. Strut collision can easily occur in tensegrity structures. Since these structures often undergo large displacement, even if the space between two struts is quite large, they may come into contact in final deformation state. There has recently been growing interest in deployable structure design based on the tensegrity concept (Korkmaz, 2011; Moored and Bart-Smith, 2009). To achieve low volume occupancy, structures may be packaged to an extreme contraction state at the folded state. In this case, it is quite reasonable to believe that one strut may contact with another one or several others. To achieve safe structural design, it is important to quantify the specific effects of collision phenomena. In tensegrity research, previous studies addressing these issues are very rare and to the best of our knowledge can be summarized as follows. Le Saux et al. (Le Saux et al., 2004) investigated the issue of a collision occurring between two slender steel bars using Moreau’s numerical model. However, their work only addressed

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

collisions occurring between two individual struts and was not extended to a whole tensegrity system. Cefalo and Mirats Tur (Cefalo and Mirats Tur, 2010) proposed two algorithms for the real-time self-collision detection of tensegrity systems by defining struts as ideal cylinders. Kan et al. (Kan et al., 2018b) found that strut collisions can frequently occur during deployment analyses of clustered tensegrity structures (Bel Hadj Ali et al., 2011; Moored and Bart-Smith, 2009; Zhang et al., 2016) when a relatively large motion speed is adopted. However, the specific effects of such phenomena have not yet been evaluated. Sultan et al. (Sultan et al., 20 0 0) proposed a flight simulator based on tensegrity structures. To enforce avoidance of strut collisions, in this work the authors maximize the distance between struts by solving a constrained quadratic optimization problem. In the robot domain, some researches (Juan and Mirats Tur, 2008; Porta and Hernández-Juan, 2016; Xu and Luo, 2012; Xu et al., 2014) have focused on the path planning of tensegrity robots with the main objective of finding a collision-free path. Other works (Juan et al., 2009) address using appropriate control techniques to achieve a collision-free path. Most of the above works simply mention collision avoidance in struts, and do not directly confront the issue of strut collision. From a practice point of view, indeed, in normal working environments strut collision are undesirable for some of these applications, as they can severely compromise system safety, potentially leading to vibrations and precision loss, component self-locking and functionality incapacitation. However, from a design point of view, it is desirable to analyze systems under various potential extreme conditions, to gain comprehensive insight into mechanical behaviors for further reliability assessment. In this sense, it is of great importance to take collisions into account in the analysis of the design phase. The classical problem of collision mechanics is quite an old topic addressed in many engineering applications and still remains as an active field of research. For modeling convenience, various approaches of analyzing the contact-collision phenomenon have already been proposed and can be generally classified into two groups: geometrical constraint-based methods (Pfeiffer and Glocker, 1996) and contact force-based methods (Lankarani and Nikravesh, 1990). For geometrical constraint-based methods, relative indentation between contacting bodies is not permitted as such bodies are considered to be entirely rigid at the contact zone. When applying these methods, a linear complementarity condition (Kwak, 1991) or differential variational inequality (Acary and Brogliato, 2008) are usually formulated based on unilateral geometrical constraints. Such methods typically require the use of well-designed algorithms (Chen et al., 2013) to obtain solutions within the mathematical framework of the measure differential inclusions, which is not an easy task. In contrast to geometrical constraint-based methods, contact force-based methods offer distinct advantages of computational simplicity and efficiency and have been extensively used in recent years (Alves et al., 2015). For such methods, two contacting bodies are permitted to penetrate slightly into one another, and the contact force present during the short-lived interval of a collision event can be computed from a continuous function of the relative penetration depth and penetration velocity. It is as if each contact region of the contacting bodies is covered with spring-damper elements (with linear or nonlinear stiffness) scattered across their surfaces. No explicit unilateral kinematic constraints are considered but simply force reaction terms are used. The normal force, including elastic and damping terms, prevents penetration. Magnitudes of the stiffness and deflection of spring-damper elements are computed based on the penetration, material properties and surface geometries of colliding bodies. The surface friction force can also be introduced using various models (Feeny et al., 1998; Oden and Martins, 1985).

203

This paper represents a first attempt at applying strut collision simulation in tensegrity mechanical analyses. The main objective is to develop a general and efficient rigid body framework based on contact force methods to analyze tensegrity systems involving strut collision. Specifically, a mathematical model of two-strut contact in three-dimensional space is addressed. A direct investigation into two representative examples is then conducted, with the main aim of revealing how collision phenomena affect structural responses and to inspire further tensegrity research by taking into account collision effects. The rest of this paper is organized as follows. Section 2 presents modeling assumptions and basic kinematic descriptions as a foundation. Then, a detailed mathematical formulation of the contact force model between two struts in threedimensional space is presented in Section 3. The related tangent stiffness and damping matrices are given in Section 4. Applications of this model to the tensegrity analyses are given in Section 5. Numerical simulations of two representative examples are presented in Section 6. Finally, conclusions are summarized in Section 7. 2. Modeling assumptions and basic kinematic descriptions Throughout this work, the following modeling assumptions are made to build a mechanical model of tensegrity structures. (1) Struts are ideal cylinders and sufficiently thick, so overall deformations can be neglected. (2) Each strut is homogeneous along its length and has a uniform cross-section. The mass distribution is uniform. (3) For class k (k > 1) (Skelton and Oliveira, 2009) tensegrity, struts are connected by frictionless ball joints. (4) The cables are very soft and light, so bending effects can be neglected. Based on the above assumptions, a strut in tensegrity systems can be considered as a rigid body in describing its configuration (Yang and Sultan, 2017). Fig. 1 shows a single rigid body located in three-dimensional Cartesian space. OXYZ denotes the global coordinate frame, and oi xi yi zi denotes the local coordinate frame fixed to the body. Generalized coordinates of the body can be defined from the configuration of the local coordinate frame:



qi = RiT

T

Θi

T

Fig. 1. A rigid body in three-dimensional Cartesian space.

(1)

204

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

where Ri ∈ R3×1 denotes the position vector of the origin of the associated local coordinate frame. Θi denotes the orientation coordinates of the local coordinate frame, which can be classically expressed by the Euler or Cardan angle (Shabana, 2005). However, these descriptive approaches may have unfavorable results such as singularities of the configuration and time-consuming calculations of trigonometric functions. To overcome these drawbacks, quaternion (Cefalo and Mirats-Tur, 2011; Shabana, 2005) is used here to describe the orientation coordinates. As a result, Θi ∈ R4×1 , and seven generalized coordinates are needed to describe the configuration of a single body (qi ∈ R7×1 ). The generalized velocities of the body can be obtained by taking the derivative of generalized coordinates with respect to time as:



q˙ i = R˙ iT

T T

Θ˙ i

(2)

An arbitrary material point P in the body can be uniquely iden− → tified by a vector u P ∈ R3×1 , the position vector of oP whose components are measured in the local coordinate frame (throughout this paper, (•) denotes a quantity (•) measured in the associated local coordinate frame). u P is a constant regardless of whether the rigid body is fixed in space or is undergoing a large overall motion. The position vector of point P measured in the global coordinate − → frame (OP ) can be calculated as:

RP = Ri + Ai u P

(3)

3×3

where Ai ∈ R is the orthogonal transformation matrix between the local coordinate frame and the global coordinate frame. The absolute velocity vector of point P can be calculated from the derivative of Eq. (3) with respect to time:

vP = R˙ P = BP q˙ s

(4)

3×7

where BP ∈ R , a widely used matrix in the following deduction, can be defined as:

BP = ∂ RP /∂ qi

(5)

It is of particular importance to know the specific structure of such a matrix. Based on the properties of orientation coordinates (Shabana, 2005), BP can be computed as



 P Li −2Ai u



3. Formulation of the contact force model In this section, a general and comprehensive mathematical model for modeling contact between two struts in threedimensional space is formulated. The model is presented first in the context of a dynamic analysis. The static analysis will be discussed as a special case. The fundamental objective is to establish the relationship between the generalized coordinates/velocities and the generalized force vector of contact forces. The detection of the collision of two struts in three-dimensional space is first addressed. 3.1. Collision detection Previous studies (Cefalo and Mirats Tur, 2010; Juan et al., 2009) detect the collision of two struts using an optimization approach. In the present work, a geometry-based method is used, and the collision condition is directly expressed using our coordinate system and variables; therefore, the results can be directly implemented in the following dynamic time history analysis or static iterative analysis. As is shown in Fig. 2, consider an imminent collision occurring between two struts. The local coordinate frame is crucial in describing the body configuration. To simplify the formulation, the origin of the local coordinate frame oi is exactly positioned at the midpoint of the corresponding axis. uˆ i , which is measured in the local coordinate frame, represents the unit vector along the axis. The length and radius of the two struts are Li and ri , respectively. Then, a natural coordinate si can be introduced to mark the arbi precisely denotes the trary material point Pi in the axes, and si u i −→ position vector of oi Pi measured in the corresponding local coordi→

nate frame. Define vector h = P2 P1 . According to Eq. (3), its components measured in the global coordinate frame can be given as:

 1 − R2 − s2 A2 u  2 h = RP1 − RP2 = R1 + s1 A1 u

(6)



 ∈ R3×3 is where I3 ∈ R is the third-order identity matrix, u P the antisymmetric matrix corresponding to vector u P , and Li ∈ 3×4 R is the transition matrix between the angular velocity and ˙ ). generalized velocity (ω i = 2Li Θ i The formulation of a concentrated force can now be presented by applying the principle of virtual work. As is shown in Fig. 1, at the current time a concentrated force is applied to the body with the action point being material point P, with the magnitude being F, and with the unit direction vector measured in the global coordinate frame being d∈ R3×1 (throughout this paper, ( · ) denotes the unit direction vector associated with vector (·)). Virtual work done by the concentrated force can be computed as



BP = I3

3×3

δW = F T δ RP

(7)

where F = F dˆ is the force vector, and δ RP is the virtual displacement of the action point P under the variation of the generalized coordinates (δ qi ). Both of these are measured in the global coordinate frame. Note from Eq. (5) that one can obtain:

δ RP = BP δ q i



 1 T h = 0 A1 u



(11)

 2 T h = 0 A2 u

Substituting Eq. (10) into Eq. (11) yields a linear equation with unknowns [s1 , s2 ]. The solutions can be computed from:

 s1 s2

=

1

 2  T AT A2 u −u 1 1

 2  T AT A2 u −u 1 1

1

−1



 T AT (R2 − R1 ) u 1 1  T AT (R1 − R2 ) u 2 2

(12)

 T AT A u  Note that Eq. (12) is invalidated when u 1 1 2 2 = ±1, which represents a configuration in which the two struts are parallel. To avoid an overly tedious discussion, in the following deductions we frequently omit such a special case. Once s1 and s2 are obtained, vector h and its norm h can be determined. h denotes the distance between the two axes of the struts. As the natural coordinate si cannot exceed the strut length, the collision condition can be finally written as: when and only when the logical relaL L L L tion(h < r1 + r2 ) ∧ (s1 ∈ [− 21 , 21 ] ) ∧ (s2 ∈ [− 22 , 22 ] ) is true. 3.2. Normal contact force vector

(8)

Substituting Eq. (8) into Eq. (7) yields the final formulation of the generalized force vector of the concentrated force:

Qc = F BTP d

(10)

√ If h = P2 P1 = hT h represents the minimum spacing between the two axes, it can be proved that the geometric perpendicular −−→ −−→ −−→ −−→ relations o1 P1 ⊥P2 P1 and o2 P2 ⊥P2 P1 must be satisfied (Schneider and Eberly, 2002). The following formula can thus be obtained. →

(9)

For contact between two general surfaces, the contact force vector can be habitually decomposed into two components: the normal force and the tangential force. The normal force prevents penetration, and the tangential force prevents relative sliding.

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

205

Fig. 2. Two struts in three-dimensional space: no collision occurs (left), collision occurs.

The contact force vector is determined by identifying its three important ingredients: magnitudes, directions and action points. The contact points of strut surfaces and the action points of the contact force, can be determined based on the previous collision −−→ detection analysis. Let h = P2 P1  represent the minimum spacing between the two strut axes. When strut collision occurs, the contact points and the material points P1 and P2 must be collinear. Let material points PA and PB be the contact points in struts 1 and 2, −→ ˆ along vector − respectively. Define a unit direction vector h P2 P1 . Relation Eq. (11) also indicates that vector  h (h) is perpendicular to the tangent plane at the contact points of the two struts ˆ denotes the direction of the normal surfaces. Therefore, vector h −→ −→ contact force. The position vectors of OPA and OPB measured in the global coordinate frame can be given as:



 1 − r1 RPA = RP1 − r1 h = R1 + s1 A1 u h



 1 − r1 h ˆ u PA = AT1 RPA − R1 = s1 u 1

(14)

 1 + r2 u PB = AT2 (RPB − R2 ) = s1 u h2

ˆ  = AT h ˆ (i = 1, 2) denotes the unit direction vector h ˆ meawhere h i i sured in the respective local coordinate frame. Using expressions ˜ ˜ Eq. (6), Eq. (14) and the tensor transformation relation hˆ = Ai ˆhi ATi , it is necessary to define several significant B matrices corresponding to the four important material points Pℵ (ℵ = A, B, 1, 2).

⎧  ⎪ ⎨BPA = I3  ⎪ ⎩BPB = I3









˜ −2A1 uˆ PA L1 = BP1 − 03×3

˜ −2r1 hˆ A1 L1

˜ −2A2 uˆ PB L2 = BP2 + 03×3

˜ −2r2 hˆ A2 L2





(15)

˜ −2si Ai uˆ i Li ] (i = 1, 2). With Eq. (4), the absolute velocity vector of material points Pℵ (ℵ = A, B, 1, 2) can be calculated as: where (using Eq. (6)) BPi = [I3



vPℵ = BPℵ q˙ 1 (ℵ = A, 1) vPℵ = BPℵ q˙ 2 (ℵ = B, 2)

(16)

At this point, the directions and the action points of the normal contact force have been fully determined, and it is now necessary to determine its magnitude, which must be computed by using a



(17)

where K, ζ , ζ˙ and D are the contact stiffness coefficient, relative penetration depth, relative collision velocity and hysteresis damping coefficient, respectively. The hysteretic damping coefficient D can be expressed as a function of the restitution coefficient ce and initial impact velocity v0 as:



D=

−−→ −−→ The position vectors of o1 PA and o2 PB measured in the respective local coordinate frames are important because they directly identify the material points at which contact takes place:





FN = K ζ 3/2 1 + Dζ˙

(13)

 2 + r2 ˆ = R2 + s2 A2 u RPB = RP2 + r2 h h



suitable constitutive law. Various types of constitutive laws have been proposed in the available literature, and the most prominent one was proposed by Hertz (Hertz, 1881). However, this law is purely elastic in nature and cannot explain energy loss occurring during the collision process. To overcome its drawbacks, Lankarani and Nikravesh (Lankarani and Nikravesh, 1990) extended Hertz’s law by taking into account the dissipative components:

3 1 − ce2



4v0

(18)

This contact force model is found to be satisfactory for general mechanical contacts, and is widely used in engineering (Tian et al., 2018). It is thus adopted here for tensegrity strut collision modeling. However, the methodology used in this work is general, and other suitable constitutive laws may also available for the normal contact force calculation. A further discussion of the differences between various constitutive laws falls beyond the scope of this paper. A recent survey of this topic can be found in (Alves et al., 2015). For contact between two struts, based on the above analysis, the relative penetration depth function can be written as:

ζ = r1 + r2 − h

(19)

It is worth emphasizing that the relative penetration is only an equivalent concept in describing the local deformation (or indentation) of two struts (Lankarani and Nikravesh, 1990; Skrinjar et al., 2018; Tian et al., 2018). The rigid body assumption is still applicable because the overall deformation is very small. The use of penetration allows one to directly measure the local deformations once the generalized coordinates are given. When strut collision occurs, ζ > 0 is always satisfied. To obtain −−→ the relative collision velocity ζ˙ , the change rate of vector h = P2 P1 can be first calculated with Eq. (16):

h˙ = vP1 − vP2 = BP1 q˙ 1 − BP2 q˙ 2

(20)

206

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

where μ is the friction coefficient and cd is the dynamic correction coefficient expressed as:

⎧ 0 ⎪ ⎪ ⎨ |vT | − v0 cd = v1 − v0 ⎪ ⎪ ⎩ 1

The length change rate of vector h, can thus been computed by ˆ: projecting h˙ to the direction vector h

(21)

Taking the derivative of Eq. (19) with respect to time and in combination with Eqs. (20) and (21), the relative collision velocity ζ˙ yields:

ζ˙ = −h˙ = −hˆ T (BP1 q˙ 1 − BP2 q˙ 2 )

(22)

Thus, the normal contact force vectors acting on the two asˆ and the corresponding sociated struts can be expressed as ±FN h generalized force vector can be computed as:



ˆ = FN BT h ˆ QN1 = FN BTPA h P1 ˆ = −FN BT h ˆ QN2 = −FN BTPB h P2

(23)

(25)

if vT > v1

vtol = vPA − vPB = BPA q˙ 1 − BPB q˙ 2

vN = hˆ hˆ T vtol

(27)

The tangential sliding velocity vector vT can be given as:



vT = vtol − vN = I3 − hˆ hˆ T vtol

(24)

(28)

The magnitude vT and the direction vˆ T can be accordingly obtained. Considering that the direction of the tangential contact/friction force vector is exactly opposite to the tangential sliding velocity, the tangential contact force vectors acting on the associated two struts can be expressed as ∓FT vˆ T . Taking into account the positions of the action points, the generalized force vector of the tangential contact force can be written in accordance with Eq. (9) as:



(29)

QT2 = FT BTPB vˆ T

The tangential contact force caused by surface friction, acts when contacting struts are prone to sliding across each other. Friction is a complex phenomenon that can lead to different friction regimes such as sticking and sliding. Generally, the most fundamental friction force model is the Coulomb’s friction law (Flores et al., 2006). This model assumes that the friction force between sliding bodies is proportional to the normal contact force. Coulomb’s friction law poses numerical difficulties when the relative tangential velocity is in the vicinity of zero. As indicated in Fig. 3(a), problems arise during integration because the friction force changes instantaneously. To overcome this drawback, in this study a modified Coulomb’s friction law (Ambrósio, 2002) is used (schematically illustrated in Fig. 3(b)). The magnitude of the tangential contact force is expressed as

(26)

Projecting the total relative velocity vector vtol to the normal ˆ yields: direction vector h

QT1 = −FT BTPA vˆ T

3.3. Tangential contact force vector

FT = μcd FN

if v0 ≤ |vT | ≤ v1

where vT is the tangential sliding velocity and v0 and v1 are the given velocity tolerances. Using Eq. (16), the total relative velocity vector observed between the material contact points PA and PB , can be given as:

Fig. 3. (a) Classical and (b) modified, versions of Coulomb’s friction law.

ˆ T h˙ h˙ = h

if vT < v0

3.4. Determination of the stiffness coefficient In this subsection, the contact stiffness coefficient involving two-cylinder contact in an arbitrary angle is given based on the Hertzian contact theory. Consider instantaneous contact illustrated in Fig. 4, where θ denotes the intersection angle and oxi yi (i = 1, 2) denotes the axis of the principal curvature of the corresponding surface. In the neighborhood of the contact point, one can found a suitable selection of axes oxy allows the separation of the two surfaces to be written as (Johnson, 1985):

h ¯ = Ax2 + By2 =

1 2 1 2 x + y 2R 2R

(30)

where A and B are positive constants. R and R  are defined as major and minor principal relative radii of curvature, respectively. The following relation holds (see details from (Johnson, 1985) in Appendix 2)

Fig. 4. The contact surfaces of two struts with its local enlarged drawing.

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

207

⎧   1 1 1 1 1 ⎪ ⎪ ⎨A + B = 2 R 1 + R 1 + R 2 + R 2   1   1   1  1  ⎪ 1 2 1 2 1 1 ⎪ ⎩B − A = 1 − + − + 2 − − cos 2θ         2

R1

R

R2

1

R

2

R1

R

R2

1

where R i and R  i (i = 1, 2) are the principal radii of the curvature of the surface of strut i. Eq. (31) is general for contact between two arbitrary surface shapes. For contact between two cylindrical struts, we have R i = ∞ and R  i = ri . Thus A and B can be determined from Eq. (31) and an equivalent radius Re can be defined as

Re =

√ R R =

1 √ 2 AB

(32)

For a given load P, the formulation of Hertzian contact theory yields (see details from (Johnson, 1985) in Appendix 3)

 δ=

9P 2 16Re E ∗2

1 / 3



F2 R /R



(33)

where δ is local deformation and E∗ is the equivalent modulus of elasticity defined as

 E∗ =

1 − v21 1 − v22 + E1 E2

−1 (34)

in which Ei and vi (i = 1, 2) are, respectively, the material elastic modulus and Poisson’s ratio associated with strut i. F2 (R /R ) is a function of the ratio of the major to minor principal relative radii of curvature. Its exact expression involves an elliptical integral and is extremely complex. As suggested in Johnson, 1985), as a first order approximation, in this work F2 (R /R ) is taken as unity. Using Eqs. (32) and ((33), the contact stiffness coefficient can be given as:

K=

P

δ

3/2

=

4E ∗  Re 3

(35)

To gain an insight into the evolution of contact stiffness, Fig. 5 plots the ratio of Kθ to K90◦ against the intersection angle θ . It can be found that for a relatively large intersection angle (θ > 60°), the contact stiffness coefficient tends to be a constant. This implies that one can use an invariable contact stiffness coefficient K90◦ for analyses of the predetermined contact of a large intersection angle, preventing the frequent updating of the stiffness coefficient.

R

(31)

2

4. Tangent stiffness matrices and tangent damping matrices The tangent stiffness matrix is defined as the derivative of the generalized force vector with respect to the generalized coordinates. The tangent damping matrix is defined as the derivative of the generalized force vector with respect to the generalized velocities. These matrices are important in ensuring convergence in static and dynamic analyses. In this section, the tangent stiffness matrix and tangent damping matrix associated with the previously derived generalized force vectors are analytically derived. 4.1. Tangent stiffness matrix of the normal contact force The tangent stiffness matrix of the normal contact force is deduced in this subsection. To simplify the derivation, this matrix is divided into four components:

⎡ ∂Q

∂ QN1 ⎤ ∂ q2 ⎥ ⎦ ∈ R14×14 ∂ QN2 ∂ q2

N1

⎢ ∂ q1 ∂ QN2 ∂ q1

KN = ⎣

(36)

Note that Eq. (15) can be first redefined for convenience,

 ⎧   ⎪ BPA = I3 03×7 +s1 03×3 ⎪ ⎪ ⎪ !" # ⎪ ⎪ ⎪ B¯ const ⎪ ⎪   ⎪ ⎪ ˜ ⎨ −r1 0 −2hˆ A1 L1 3×3    ⎪ I 0 B = + s ⎪ 3 3 ×7 P 2 03×3 B ⎪ ⎪ !" # ⎪ ⎪ ⎪ B¯ const ⎪ ⎪   ⎪ ⎪ ⎩ +r2 0 ˜ −2hˆ A L 3×3

˜ −2A1 uˆ 1 L1

!"



#

B¯ 1

˜ −2A2 uˆ 2 L2

!"



(37)

#

B¯ 2

2 2

and hence

BPi = B¯ const + si B¯ i (i = 1, 2)

(38)

It should be very clear from Section 3.1 (Eq. (12)) that either of the two natural coordinates (s1 /s2 ) is a function of both q1 and q2 . Taking the partial derivatives of Eq. (10) with respect to the generalized coordinates, we have

⎧ ∂h ∂ s1 ∂ s2 ⎪ = BP1 + A1 uˆ 1 − A2 uˆ 2 ⎨ ∂ q1 ∂ q1 ∂ q1 ⎪ ⎩ ∂ h = A1 uˆ 1 ∂ s1 − BP − A2 uˆ 2 ∂ s2 2 ∂ q2 ∂ q2 ∂ q2

(39)

Taking the partial derivatives of perpendicular Eq. (11) with respect to the generalized coordinates and in combination with Eq. (39) we have a linear algebraic equations of ∂ si /∂ qj (i, j = 1, 2) whose solution can be written as:

⎡ ∂s

∂ s1 ⎤ 1 ∂ q2 ⎥ ⎦= T ∂ s2 −uˆ 1 AT1 A2 uˆ 2 ∂ q2

1

⎢ ∂ q1 ⎣ ∂ s2 ∂ q1

T

−hT B¯ 1 − uˆ 1 AT1 BP1 T

uˆ 2 AT2 BP1 Fig. 5. The evolution of the ratio of Kθ to K90◦ against the intersection angle.

T

uˆ 1 AT1 BP2 T

hT B¯ 2 − uˆ 2 AT2 BP2

T

−uˆ 1 AT1 A2 uˆ 2

−1

1

 (40)

208

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

Then Eq. (39) is completely determined, and some intermediate terms can then be derived from Eqs. (19) and (22).

√

hT ∂ h ∂ζ ∂h ∂ hT h =− =− =− ∂ q1 ∂ q1 ∂ q1 h ∂ q1

(41)

(42)







∂ q˙ T1 BTP1 hˆ − q˙ T2 BTP2 hˆ ∂ BTP1 hˆ ∂ BTP2 hˆ ∂ ζ˙ =− = −q˙ T1 + q˙ T2 ∂ q1 ∂ q1 ∂ q1 ∂ q1

(43)

Taking the derivative of Eq. (17) with respect to the generalized coordinates of strut I and using Eqs. (41) and (43) yields:



∂ FN ∂ζ = K 1 + Dζ˙ ζ 1/2 + KDζ 3/2 ∂ q1 2 ∂ q1 ∂ q1

∂ ζ˙

(44)

A similar expression for ∂ FN /∂ q2 can also be obtained. Here we omit variations caused by the intersection angle θ and let ∂ K /∂θ = 0 → ∂ K /∂ qi = 0(i = 1, 2), as it is reasonable to believe that the stiffness coefficient does not change considerably in each iteration step of analyses. This is even more precise in the case of contact with a large intersection angle as illustrated in Section 3.4. Next, noting from Eq. (37) that B¯ 1 /B¯ 2 is a matrix independent of q2 /q1 and using Eq. (38) we have

⎧ T

∂ BP1 hˆ ⎪ ⎪ ⎪ ⎪ ⎪ ∂ q1 ⎪ ⎪ ⎪ ⎪

⎪ ⎪ ⎪ ∂ BTP1 hˆ ⎪ ⎪ ⎨ ∂ q2 T

⎪ ∂ BP2 hˆ ⎪ ⎪ ⎪ ⎪ ⎪ ∂ q1 ⎪ ⎪ T

⎪ ⎪ ˆ ⎪ ⎪ ⎪ ∂ BP2 h ⎪ ⎩ ∂ q2

=

BTP1



$ ∂ B¯ T1V3 $$ ∂ hˆ + s1 ∂ q1 ∂ q1 $$

ˆ + B¯ T1 h

ˆ V3 =h

= BTP1

∂ s1 ∂ q1

∂ hˆ ˆ ∂ s1 + B¯ T1 h ∂ q2 ∂ q2

∂ hˆ ˆ ∂ s2 + B¯ T2 h ∂ q1 ∂ q1 T $$ ∂ B¯ 2V3 $ ∂ hˆ = BTP2 + s2 $ ∂ q2 ∂ q2 $

ˆ V3 =h

ˆ + B¯ T2 h

∂ s2 ∂ q2



∂ BTP1 hˆ ˆ ∂ FN = FN + BTP1 h ∂ q1 ∂ q1 T

∂ BP1 hˆ ˆ ∂ FN = FN + BTP1 h ∂ q2 ∂ q2

∂ BTP2 hˆ ˆ ∂ FN = −FN − BTP2 h ∂ q1 ∂ q1

∂ BTP2 hˆ ˆ ∂ FN = −FN − BTP2 h ∂ q2 ∂ q2

⎢ ∂ q˙ 1 ∂ QN2 ∂ q˙ 1

DN = ⎣

∂ QN1 ⎤ ∂ q˙ 2 ⎥ ⎦ ∈ R14×14 ∂ QN2 ∂ q˙ 2

(47)

Based on Eqs. (17) and (22), we have

⎧ ∂ FN ∂ ζ˙ ⎪ ⎪ ⎨ ˙ = KDζ 3/2 ˙ = −KDζ 3/2 hˆ T BP1 ∂ q1 ∂ q1 ˙ ⎪ ⎪ ⎩ ∂ FN = KDζ 3/2 ∂ ζ = KDζ 3/2 hˆ T BP2 ∂ q˙ 2 ∂ q˙ 2

(48)

Using Eq. (23), the tangent damping matrix of the normal contact force can be finally written as:

⎧ ∂Q N1 ⎪ ⎪ ⎪ ∂ q˙ 1 ⎪ ⎪ ⎪ ⎪ ∂ QN1 ⎪ ⎪ ⎨ ˙ ∂ q2 ∂ QN2 ⎪ ⎪ ⎪ ⎪ ∂ q˙ 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ QN2 ∂ q˙ 2

∂ FN ˆh ˆ T BP = −KDζ 3/2 BTP1 h 1 ∂ q˙ 1 ˆ ∂ FN = KDζ 3/2 BT h ˆ ˆT = BTP1 h P1 h BP2 ∂ q˙ 2 ˆ ∂ FN = KDζ 3/2 BT h ˆ ˆT = −BTP2 h P2 h BP1 ∂ q˙ 1 ˆ ∂ FN = −KDζ 3/2 BT h ˆ ˆT = −BTP2 h P2 h BP2 ∂ q˙ 2 ˆ = BTP1 h

(49)

4.3. Tangent stiffness matrix of the tangential contact force

= BTP2



N1

(45)

In accordance with Eq. (23), the block matrix in Eq. (36) can be written as:

⎧ ⎪ ∂ QN1 ⎪ ⎪ ⎪ ⎪ ∂ q1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ∂ QN1 ⎪ ⎨ ∂ q2 ⎪ ⎪ ∂ QN2 ⎪ ⎪ ⎪ ⎪ ∂ q1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ QN2 ∂ q2

Since dissipation terms are introduced in the contact force model, the tangent damping matrix must also be evaluated. Similarly, this matrix is divided into four components:

⎡ ∂Q

  hhT ∂ h ∂ hˆ ∂ (h/h ) 1 = = I3 − 2 ∂ q1 ∂ q1 h ∂ q1 h

3

4.2. Tangent damping matrix of the normal contact force

The tangent stiffness matrix of the tangential contact force is given in this subsection. This matrix, however, is more complex than the other matrices derived above. To avoid overly cumbersome operations, some necessary approximations are introduced. To simplify the derivation, the matrix can also be divided into four components similar to those of Eq. (36). Some intermediary and key terms can be first obtained. Similar to Eq. (42), we have

  vT vT ∂ vT ∂ vˆ T ∂ (vT /vT ) 1 = = I3 − 2 T ∂ q1 ∂ q1 vT vT ∂ q1

(50)

ˆ T vtol = −ζ˙ . Using Eqs. (22) and (26), it can be verified that h Thus,

(46)

By using Eqs. (44) and (45), the final expression of the tangent stiffness matrix can be computed. It can be observed that to obtain the exact tangent stiffness matrix, one should also evaluate terms such as ∂ (B¯ T1V3 )/∂ q1 , etc. Theoretically, these terms can be written in the form of a third-order tensor by multiplying a corresponding vector. However, this operation is complex and not conducive for designing programs. To simplify numerical calculations, the detailed expressions of these terms are directly given in a component manner in the Appendix A of this paper.



˙ ∂ vT ∂ vtol ∂ hˆ ζ˙ ∂ vtol ˙ ∂ hˆ ˆ ∂ζ = + = +ζ +h ∂ q1 ∂ q1 ∂ q1 ∂ q1 ∂ q1 ∂ q1

(51)





∂ B¯ 1 q˙ 1 ∂ vtol ∂ BPA q˙ 1 − BPB q˙ 2 ∂ s1 ∂ s2 ¯ = ≈ B1 q˙ 1 + s1 − B¯ 2 q˙ 2 ∂ q1 ∂ q1 ∂ q1 ∂ q1 ∂ q1 (52) Substituting Eq. (51) into Eq. (50) yields the final expression for

∂ vˆ T /∂ q1 . A similar expression can be obtained for ∂ vˆ T /∂ q2 . Note from Eq. (24) that we have

∂ FT ∂ FN = μc d ( i = 1, 2 ) ∂ qi ∂ qi

(53)

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

In accordance with Eq. (29), the tangent stiffness matrix of the tangential contact force can be computed from

⎧ ⎪ ⎪ ∂ QT1 ⎪ ⎪ ∂ q1 ⎪ ⎪ ⎪ ⎪ ⎪ ∂ QT1 ⎪ ⎪ ⎨ ∂ q2 ⎪ ∂ QT2 ⎪ ⎪ ⎪ ⎪ ⎪ ∂ q1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ QT2 ∂ q2



∂ BTPA vˆ T ∂ FT = −FT − BTPA vˆ T ∂ q1 ∂ q1

∂ BTPA vˆ T ∂ FT = −FT − BTPA vˆ T ∂ q2 ∂ q2

∂ BTPB vˆ T ∂ FT = FT + BTPB vˆ T ∂ q1 ∂ q1

∂ BTPB vˆ T ∂ FT = FT + BTPB vˆ T ∂ q2 ∂ q2



qTN

T

(59)

Since quaternion is used as the orientation coordinates, there are 7N generalized coordinates altogether. In accordance with basic theories of mechanical mechanics such as the Lagrangian equation of motion the governing equation of the system can be derived as the following index-3 differential algebraic equations (DAEs) (Shabana, 2005):

M (q )q¨ +Φq λ+Qi (q, q˙ ) − Qe (q, q˙ , t ) = 0

⎧

T $$ T ˆ ∂ B v B¯ 1V3 $ ∂ ⎪ ˆ T ∂ v ∂ s1 PA T ⎪ ⎪ ≈ BTPA + s1 + B¯ T1 vˆ T $ ⎪ ⎪ ∂ q ∂ q ∂ q ∂ q1 $ 1 1 1 ⎪ ⎪ V3 =vˆ T ⎪ ⎪ ⎪ ⎪ ∂ BT vˆ

⎪ ∂ vˆ T ¯ T ∂ s1 ⎪ PA T ⎪ ≈ BTPA + B1 vˆ T ⎨ ∂ q2 ∂ q2 ∂ q2 T

⎪ ˆ ∂ BPB vT ∂ vˆ T ¯ T ∂ s2 ⎪ ⎪ ≈ BTPB + B2 vˆ T ⎪ ⎪ ∂ q ∂ q1 ∂ q1 ⎪ 1 ⎪ ⎪

T $$ ⎪ ⎪ T ⎪ ∂ s2 ⎪ ∂ BB vˆ T ≈ BT ∂ vˆ T + s ∂ B¯ 2V3 $ ⎪ + B¯ T2 vˆ T ⎪ 2 PB ⎩ ∂ q2 ∂ q2 ∂ q2 $$ ∂ q2

(55)

4.4. Tangent damping matrix of the tangential contact force The tangent damping matrix of the tangential contact force is given in this subsection. Similar to Eq. (47), this matrix is divided into four components for convenience. One must first note from Eqs. (24) and (48) that

⎧ ∂ FT ∂ FN ⎪ ⎨ ˙ = μcd ˙ = −μcd KDζ 3/2 hˆ T BP1 ∂ q1 ∂ q1 ⎪ ⎩ ∂ FT = μcd ∂ FN = μcd KDζ 3/2 hˆ T BP 2 ∂ q˙ 2 ∂ q˙ 2

T

Φ ( q, t ) = 0

V3 =vˆ T

(56)

From Eqs. (26) and (28), it can be verified that

   

vT vTT ∂ vT vT vTT ∂ vˆ T 1 1 ˆ T BP (57) ˆh = I − 2 = I3 − 2 I3 − h A ∂ q˙ 1 vT 3 vT ∂ q˙ 1 vT vT A similar expression can be obtained for ∂ vˆ T /∂ q˙ 2 . With Eq. (29) the tangent damping matrix of the tangential contact force can be computed from

∂ vˆ T ∂ FT − BTPA vˆ T ∂ q˙ 1 ∂ q˙ 1 ∂ vˆ T ∂ FT = −FT BTPA − BTPA vˆ T ∂ q˙ 2 ∂ q˙ 2 ∂ vˆ T ∂ FT = FT BTPB + BTPB vˆ T ∂ q˙ 1 ∂ q˙ 1 ∂ vˆ T ∂ FT = FT BTPB + BTPB vˆ T ∂ q˙ 2 ∂ q˙ 2

···

qT2



where

⎧ ∂ QT1 ⎪ ⎪ ⎪ ⎪ ∂ q˙ 1 ⎪ ⎪ ⎪ ⎪ ∂ QT1 ⎪ ⎪ ⎨ ˙ ∂ q2 ⎪ ∂ QT2 ⎪ ⎪ ⎪ ⎪ ∂ q˙ 1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ QT2 ∂ q˙ 2

rigid bodies. For a structure containing N struts, the generalized coordinates of the system are directly chosen as an assemblage of generalized coordinates of each strut:

q = qT1

(54)

209

(60)

where M is the mass matrix, Φ is the constraint equation and Φq is the Jacobian matrix of the constraint equations, λ is the Lagrange multiplier vector physically representing the constraint forces, qi is the quadratic velocity vector (the sum of Coriolis and centrifugal force). Expressions of the mass matrix M and the quadratic velocity vector qi can be found in (Shabana, 2005). qe is the generalized external force vector, which is induced by various kinds of applied force such as typical external concentrated forces, elastic forces of classical and clustered (Bel Hadj Ali et al., 2011; Moored and Bart-Smith, 2009) cables and, the main concerns of this work, strut collision forces. For the latter ones, relationships between the generalized coordinates/velocities and generalized force vectors are given in Section 3. For a typical tensegrity system, constraint equation Φ may include three parts. The first one is induced by the orientation coordinates of each strut due to the use of quaternion written as

% %2 Φ ( 1 ) = %Θ i % − 1 = 0

(61)

where  ·  denotes the 2-norm of the associated vector. The related constraint Jacobian matrix can be given (at a single strut level) as

∂ Φ (1 )  = 0 ∂ qi

0

0

2Θ i

T



(62)

The second one is induced by joints between different struts. For a class k (k ≥ 2) tensegrity system, two or more struts may be connected at their ends by ball joints as shown in Fig. 6. From the kinematic descriptions presented in Section 2, the constraint equations can be easily expressed as:

Φ ( 2 ) = R1 + A1 u 1 − R2 − A2 u 2 = 0

(63)

Accordingly, the associated constraint Jacobian matrix at the joint element level can be given as

= −FT BTPA

Φq(2e ) = (58)

$  ∂ Φ (2 )  = I3 −2A1 u˜ 1 L1 $−I3 2A2 u˜ 2 L2 ∂ qe

where qe = [qT1 qT2 ]T is the assemblage of the generalized coordinates of the related struts.

Clearly, all matrices derived in this section are only applicable for a case in which a collision occurs. Otherwise, these terms should been set as zero matrices. 5. System modeling, dynamic equation and numerical solving Based on the assumptions described above, a tensegrity system can be modeled as a general multibody system in which the struts are considered as rigid bodies, and the elastic internal force of the cables can be considered as a set of applied forces acting on the

(64)

Fig. 6. Two struts connected with a ball joint.

210

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

The previous two parts are holonomic and scleronomic constraints while the third part represents a set of rheonomic constraints, i.e., it has explicit time dependence. As mentioned previously, there has been growing interest in deployable structure design using the tensegrity concept. To deploy or fold structures, some components of the system may be subjected to displacement motions. This part of constraint equations is specially introduced to simulate the actuation behavior of an active tensegrity system, which can be written in a general mathematical form

Φ ( 3 ) = g ( q, t )

(65)

Detailed expressions of the constraint equations and the associated constraint Jacobian matrix, depend on artificial settings or specific control algorithms. With regard to the numerical solving of dynamic equations, many direct integration methods have been proposed to solve index-3 DAEs. In this study, we use the classical Newmark method (Newmark, 1959), which assumes that:



q˙ n+1 = q˙ n + [(1 − δ )q¨ n + δ q¨ n+1 ]t

(66)

qn+1 = qn + q˙ n t + [(1 − 2α )q¨ n + 2α q¨ n+1 ]t 2 /2

where t is the integration step size, (·)n and (·)n + 1 are the terms corresponding to the nth and n + 1th time steps, respectively, and α and δ are the Newmark parameters. Discretized equilibrium equations at the n + 1th time step yield

fn+1 =



Mn+1 q¨ n+1 + Φqn+1 λn+1 + Qn+1 (qn+1 , q˙ n+1 , tn+1 ) T

Φ(qn+1 , tn+1 )/α t 2

=0 (67)

where the term q denotes the total generalized force (qi − qe ) in Eq. (60). The constraint equations Φ multiplied by a scaling factor 1/α t2 increase the stability of the algorithm (Bottasso et al., 2007). While adopting X = [q¨ Tn+1 λTn+1 ]T as unknowns in combination with Eq. (66), the solution of Eq. (67) can be obtained using the Newton–Raphson method with the Jacobian matrix for each time step given as

⎡

J=

∂ fn+1 ⎢Mn+1 + ∂ (Mn+1 q¨ n+1 ) + α t 2 =⎣ ∂ qn+1 ∂X Φqn+1

  ∂ ΦTq λn+1 ∂ qn+1

+ α t 2

It can be clearly observed that to obtain the exact Jacobian matrix one must also evaluate terms such as ∂ q/∂ q and ∂ Q /∂ q˙ . For contributions induced by strut collision forces, these terms are the tangent matrices derived in Section 4. For contributions induced by cable forces, these terms can be found in (Kan et al., 2018b). Fig. 7 shows a computational chart for numerically solving the governing equations. Note that cable slacking can occur in the analysis. In this work, cable slacking is handed in a straightforward manner: if it is found that in the current iteration some cables are slacking, in the next iteration these cables will not contribute to the system generalized force vector and the Jacobian matrix. The above sections focus on analyses conducted from a dynamic point of view. Extensions to static analysis are straightforward. This can be done with two approaches. The first one is conceptually similar to the conventional dynamic relaxation method (Barnes, 1999), i.e., by introducing pseudo damping terms or other dissipation forces and by performing a conventional dynamic analysis; the system will eventually evolve into a static balance configuration. The second one is a direct method: by directly letting the generalized velocity and the generalized acceleration in Eq. (60) to be zero, the dynamic governing equation can degenerate into a set of nonlinear algebraic equations representing the static equilibrium

Table 1 The material and dimension parameters. Parameter

Value

Section radius of the cables Young’s modulus of the cables Length of the struts Section radius of the struts Young’s modulus of the struts Poisson’s ratio of the struts Prestress force in the lower cables Prestress force in the upper cables Prestress force in the bracing cables Prestress force in the struts

1 cm 7.84 MPa 122.47 cm Rs 210 GPa 0.31 200 N 282.84 N 282.84 N 489.89 N

conditions. In this case, the tangential contact force vector cannot be modeled due to the use of the modified Coulomb’s friction law (Eq. (24)), and the formulation of the normal contact force degenerates into Hertz’s law according to Eq. (17). Solutions of the static analysis can be obtained by using the Newton–Raphson method with an appropriate selection of initially estimated solutions. This approach is adopted in our following static analyses. 6. Numerical simulations In this section, two representative numerical examples involving strut collision phenomena are presented. The purpose is, for one hand, to demonstrate the effectiveness of the general framework, in particular the contact force model; for the other hand to show how strut collisions influence the mechanical responses of the systems, which further enlarges the known vision of the behaviors of such structures. 6.1. A quadruplex tensegrity unit The first example involves a very classical quadruplex tensegrity unit that has been broadly used as an elementary cell in designing large-scale grid structures (Kan et al., 2018a; Zhang et al., 2015).

⎤ ∂ Qn+1 ∂ Qn+1 + δ t ∂ qn+1 ∂ q˙ n+1

ΦTqn+1 ⎥ ⎦

(68)

0 While the mechanical properties of this unit have been extensively investigated, to the best of our knowledge, no studies have addressed strut collision effects. The main focus of this study is to test the mechanical properties of such a unit under quasi-static axially compressed loading. Fig. 8 presents a perspective and a top view of the unit, where grey thick cylinders denote the struts and thin pink lines denote cables. More information is given in Table 1. The four lower points P1-P4 are fixed as boundary conditions. The four upper points P5-P8 are the loading points; each one is subjected to actuation to achieve a specified displacement control. The actuations are unified so that the surface enclosed by the four upper cables remains parallel to the x–y plane. The actuation value υ is defined as the vertical downward moving distance. It is also assumed that the actuation speeds are extremely slow, so inertia effects can be neglected. First, as a preliminary test to analyze the general structural behavior, strut collisions are not considered, as is often the case in the existing literature. The system thus evolves as if the struts can penetrate each other. The actuation is deliberately performed until a large actuation value is achieved. Fig. 9 presents the resultant of the four actuation forces vs. the actuation value curve,

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

211

Fig. 7. Computational chart for numerically solving the governing equations.

Fig. 8. The quadruplex tensegrity unit: a perspective view (left) and a top view.

which clearly exhibits axial stiffness. For a broad range of axial loading values, the structure gradually exhibits a stiffeningsoftening-stiffening behavior, which agrees with outcomes given in (Fraternali et al., 2015). Numerical results also show that an ultimate configuration exists, which corresponds to an actuation value of υ = 1.72 m. In such a configuration, actuation forces tend to infinity and the four struts are almost vertical. This is reasonable due to the rigid body assumptions of the struts. Fig. 9 also presents two typical configurations corresponding to actuation values of 0.5 m and 1 m. It is clear that strut collisions take place very early on: no later than the actuation value of 0.5 m, where the structure is compressed to a height of zero. Apparently, all subsequent simulation results are not real in practice. In what follows, we only focus on actuation value interval υ ∈ [0, 0.5]. Next, strut collisions are taken into account in our analyses. Analyses are carried out for five cases: the section radius of struts

Rs ranges from 4 to 8 cm with a step size of 1 cm. Fig. 10 gives the resultant of the actuation forces vs. the actuation value curves. It can be observed that a rich variety of behaviors emerges when considering strut collision effects. In the first stage, it is as expected that the curves merely follow the results of the previous case involving no collision. Since the actuation value is small, gaps between the struts are large enough and collisions do not really occur. As the actuation value increases, such gaps decrease gradually (as is clearly indicated in Fig. 11), and the struts begin to touch. This results in a dramatic increase in the actuation force and in axial stiffness. Before and after the collisions occurring, the structure is characterized by a very soft and a very stiff deformation mode, respectively. The section radius of the struts directly determines when collision occurs: the smaller the radius, the later the collision occurs. Once two struts touch, they slide obliquely downward along the surfaces. One weakness of the present work is that,

212

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

Fig. 9. The resultant of the actuation force vs. the actuation value.

for the static analysis, the sliding friction effect cannot be modeled, due to the modified Coulomb’s law used in this work. It can also be observed from Fig. 10 that an ultimate axial stiffness exists for each case, and interestingly the maximum value is closely related to the section radius of the struts. This ultimate state corresponds to a configuration in which a strut is just about to slip out of the length range of the other strut. To verify the modeling assumption, further calculations show that the maximum bending and axial deformation of the most slender case (Rs = 4 cm) are, respectively, of 1 × 10−4 and 1 × 10−7 m order of magnitude, demonstrating the applicability of the present rigid body assumption. Some snapshots of representative configurations Rs = 6 cm are shown in Fig. 12, where green lines denote slack cables. Cable slacking is not observed in the previous non-collision-considered case. This simple example highlights that when taking into account strut collision effects, structural behaviors can be significantly changed. The outcomes presented here thus provide important inspiration for further structural designs using such kind of tensegrity unit, by fully taking advantage of collision effects. For example, it is conceivable that by deliberately design strut geometry parameters, structural deformation modes (extreme hard/soft behaviors) can be precisely regulated, which can be further exploited to develop special mechanical materials such as extremal materials (Milton and Cherkaev, 1995) or metamaterials (Al SabouniZawadzka and Gilewski, 2018) design. 6.2. A four-stage tensegrity tower

Fig. 10. The actuation force vs. the actuation value.

Fig. 11. Locally enlarged drawing of the actuation force vs. the actuation value.

A more complex example of folding analyses of a clustered tensegrity tower is given in this subsection. Similar structures have been studied in (Bel Hadj Ali et al., 2011; Kan et al., 2018b; Zhang et al., 2016). These works are conducted, however, mostly from a quasi-static perspective, except for the one presented in (Kan et al., 2018b), for which the dynamic deployment performance is investigated. It has been demonstrated that strut collisions can frequently occur when an oversized actuation speed is applied. However, the specific effects of such phenomena have not yet been evaluated due to a lack of effective analytical tools available. In the present work, both quasi-static deployment analysis and dynamic deployment analysis are carried out to thoroughly examine effects of strut collisions. The tower is composed of an assembly of four identical quadruplex units configured with twelve strut-to-strut connections, rendering it a class 2 tensegrity structure. A detailed description of this construction can be found in (Kan et al., 2018b). Fig. 13 presents three views of the structure, where thick grey cylinders denote struts, thin pink lines denote classical cables and blue lines denote the sliding cables. The sliding cables are connected to four additional bodies (B1–B4). Thus, actuations of the external bodies are completely equivalent to actuations of the sliding cables. The tensegrity tower has a width/length of 1 m and a height of 2 m in the initial fully unfolding state. All struts are made from steel with surfaces covered with slightly thicker rubber, which is specially introduced to extend buffering time and reduce contact stiffness. All cables are made from fibreglass. The middle eight short classical cables are replaced with light spring elements with a constant axial stiffness of 20 0 0 N/m. Instead of using original classical cables with an axial stiffness of the level of 1 × 106 N/m, this modification is introduced to fold/unfold the tower using less energy, as these elements serve as the main deformation components at the folding state. More information on the structure is given in Table 2. At the initial state, no external or prestress force is applied to the structure, and all components are in their free length state. The bottom four corner points (P1–P4) are fixed on the ‘ground’. Then, cable-based actuations are conducted to fold the structure

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

213

Fig. 12. Snapshots of configurations corresponding to actuation values of (a) 0.1, (b) 0.2, (c) 0.3, and (d) 0.4 m.

Fig. 13. (a) Perspective, (b) side and (c) top views of the structure.

in such a way that the sliding cables are simultaneously pulled by driving the four additional bodies (B1–B4) at a unified speed V until a designated actuation value is achieved. The actuation value υ is defined as the moving distance of the additional bodies. For a deployment process, a reverse procedure can be performed. In this study, only the folding process is investigated. The objective is to fold such a structure to a height of less than one-quarter of its normal height, 0.5 m, to meet the specific volume requirement. To

simulate the zero-gravity environment of outer space, self-weight is not taken into account in the analyses. 6.2.1. Quasi-static folding analysis First, a quasi-static folding analysis is conducted with the aim of identifying basic characteristics of such a structure, and to help determine an appropriate actuation value to achieve the height requirement. It is assumed that the driving speeds of the four ad-

214

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219 Table 2 Parameters associated with the structure. Parameter

Value

Section radius of the cables Young’s modulus of the cables Diameter of the outer layer of the struts Diameter of the inner layer of the struts Poisson’s ratio of the outer layer of the struts Young’s modulus of the outer layer of the struts Young’s modulus of the inner layer of the struts Density of the cables Density of the outer layer of the struts Density of the inner layer of the struts Velocity tolerances v0 Velocity tolerances v1 Friction coefficient Restitution coefficient

3 mm 73 GPa 8 cm 2.5 cm 0.47 7.84 MPa 210 GPa 2450 kg/m3 930 kg/m3 7800 kg/m3 0.001 m/s 0.05 m/s 0.1 0.5

ditional bodies (B1–B4) are extremely slow, so inertia effects can be neglected. Solutions can thus be obtained by directly solving the static equilibrium equations. Note that theoretically, C216 = 120 possible contact pairs exist. However, according to a preliminary analysis, one can reduce the number of possible contact pairs to sixteen, that is, only the contacts pairs involving struts in the same layer are considered. Fig. 14 illustrates the evolution of structural height with respect to the actuation value. As a comparison, a case wherein strut collisions effects are not taken into account is also analyzed. It can be observed that as the actuation value increases, the structural height decreases in both cases, demonstrating that the tower is being folded gradually. In the initial folding phase, the two evolution curves are completely overlapped, indicating that strut collisions will not occur. Once the actuation value extends beyond this scope, considerable differences emerge: for the no-collisionconsidered case, the structural height rapidly becomes zero within a very limited actuation value increment and remains constant (zero) thereafter, while for the collision-considered case, the structural height decreases slowly with an increasing actuation value. It is clear that the no-collision-considered case yields an unreal result. The critical point, υ = 0.638, exactly denotes a configuration wherein the struts are just about to collide, which corresponds a structural height of 0.799 m. Since the objective is to fold such a tower to a height of less than 0.5 m, actuation should be continued even as the struts begin to collide. As indicated in Fig. 14, the

Fig. 14. The evolution of structural height with respect to the actuation value.

Fig. 15. Time history curves of the heights of the four layers for V = 5 m/s.

minimum structural height requirement corresponds to an actuation value of υ = 1.448. To allow for some redundancies, an actuation value of υ = 1.5 is chosen as the designated maximum actuation value for the following dynamic analyses; the configuration corresponding to such an actuation value is set as the designed final folding state. For such a state, calculations show that the maximum bending and axial deformation of the struts are valued at 1 × 10−6 m orders of magnitude, demonstrating the applicability of the present rigid body assumption. Further calculations also show that all of the components work within the permissible range at the folding state, which shows that even when strut collision occurs, structural safety can still be ensured in this particular case. 6.2.2. Dynamic folding analysis The above quasi-static analysis illustrates general features of structure folding. However, physically quasi-static actuation is only an abstract concept. It is unclear to what extent actuation can be regarded as quasi-static. To satisfy quasi-static conditions as much as possible, the actuation speed must be very low, which results in an extremely long actuation time that may compromise maneuverability. In this paragraph a dynamic folding analysis is carried out to reveal differences with the quasi-static results and identify a limit actuation speed. Dynamic actuation is applied in a manner similar to that of the quasi-static actuation but with a certain speed of V. Five actuation speeds covering different levels are considered: V = 5, 2, 0.5, 0.05, and 0.005 m/s. Since the total actuation value is 1.5 m, corresponding folding times are 0.3, 0.75, 3, 30, and 300 s, respectively. Dynamic solutions are obtained by the Newmark method using parameters α = 0.6, δ = 0.3. For all the analyses, the time step sizes and convergence errors are taken as 2 × 10−4 and 1 × 10−6 s, respectively. The analysis time for the first two cases is set to 6 s, while that for the last three cases are set to 8, 40, and 330 s, respectively, given their long actuation time. The analysis time extending beyond the actuation time is used to observe structural vibration characteristics. Figs. 15–19 present time history curves of the heights of the four layers for each case. Corresponding quasi-static results are also presented for comparison. It is clear that the height curves vary. For a very low actuation speed (V = 0.005 m/s), the height curves of all layers are extremely close to the quasi-static solutions for all analysis time. In initial phases of actuation, the height curves of all layers decrease with an equal proportion

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

Fig. 16. Time history curves of the heights of the four layers for V = 2 m/s.

Fig. 17. Time history curves of the heights of the four layers for V = 0.5 m/s.

Fig. 18. Time history curves of the heights of the four layers for V = 0.05 m/s.

215

Fig. 19. Time history curves of the heights of the four layers for V = 0.005 m/s.

of the layer index, indicating that the structures is being folded with uniform contraction. Once the designated actuation value is achieved, the height curves of all layers slightly vibrate along their target heights. For a slightly higher actuation speed (V = 0.05 m/s), height curves generally follow the quasi-static solutions. However, differences between the dynamics and the quasi-static solutions are much more notable in the time interval prior to the first collision. Once struts begin to touch, the differences soon decline to an inconspicuous level. This leads us to an important conclusion: strut collision behavior tends to decrease the amplitude of structural vibrations. Therefore, allowing for strut collisions in folded states may enhance structure shock-absorbing. For a higher actuation speed (V = 2, and 0.5 m/s), although the structure will eventually be folded to the target height, motion behaviors are much more complex during the intermediate process. In the initial phase of actuation, motion paths travel straight down due to rapid actuation behavior. Then, the height curves of all layers stage several strong rallies caused by strut collisions. For a very high actuation speed (V = 5 m/s) the overall structural height (the 4th layer) decreases rapidly with an increasing actuation value, but the height of the lower (1st and 2nd) layers experience some slightly increasing at first. This results in an extremely non-uniform contraction of the structure. Fig. 20 presents some snapshots that clearly demonstrate this point. It is reasonable to believe that in this case dynamic effects are sufficiently notable. To clearly demonstrate collisions effects, Figs. 21–25 present the time history curves of the normal collision force for the five cases. Note that due to symmetries, in the dynamics analysis the collision forces of contact pairs located in the same layer are exactly the same, while for the quasi-static analysis the collision forces of all contact pairs are the same. It can be clearly observed that the actuation speed has a strong influence on the strut collision force: a high actuation speed can give rise to a high force peak value during the folding process; the magnitude is far greater than that of the finally folded state. With a relatively low actuation speed (V = 0.05 and 0.005 m/s), the collision forces of different layers tend to be consistent and are fairly close to the quasi-static analysis results. Slight vibrations emerge at two time points: the first one corresponds to an event during which struts begin to touch, and the second one corresponds to an event during which actuation begin to terminate. Vibrations may be further avoided by applying a smooth actuation speed throughout the whole process (i.e., no sudden stopping occurs). From the above figures it seems that

216

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

Fig. 20. Snapshots of the configurations for case V = 5 m/s.

Fig. 21. Time history curves of the normal collision force for case V = 5 m/s.

Fig. 22. Time history curves of the normal collision force for case V = 2 m/s.

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

Fig. 23. Time history curves of the normal collision force for case V = 0.5 m/s.

Fig. 26. Comparison of the maximum tangential contact force values.

Fig. 24. Time history curves of the normal collision force for case V = 0.05 m/s.

Fig. 27. Comparison of maximum stress values of struts for all cases.

Fig. 25. Time history curves of the normal collision force for case V = 0.005 m/s.

217

struts in the first (bottom) layer are generally subjected a stronger impact in most cases relative to struts of the other (upper) layers. Fig. 26 compares maximum tangential contact forces for the five cases. It can be observed that the tangential forces are also closely related to the actuation speed. For a particular actuation speed, the maximum tangential force is far lower than the maximum normal force, indicating that total collision forces are dominated by the normal force. In addition, it is important to conclude that in contrast to height (position) curves, both kinds of collision forces are very sensitive to the actuation speed. Finally, our goal is to determine a level of the limit actuation speed purely from a material failure point of view. We first calculate the maximum bending stress of the struts for all cases, as shown in Fig. 27. It can be clearly observed that the maximum stress values exceed the material ultimate stress for the first two cases, which indicates that the struts are likely to break in these cases. Thus, the limit actuation speed may be valued at V = 0.5 m/s. With this actuation speed we then calculate the maximum compressive force of the struts and the maximum stress values of the classical and the sliding cables. These values are much lower than the respective allowable limits, indicating that all of the com-

218

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

ponents work within the permissible range. The above facts finally demonstrate that the limit actuation speed lie in the level of V = 0.5 m/s. In the quickest way, the system can be actuated in only several seconds; thus, a rapid actuation behavior can still be achieved.

Dr. Lingchong Gao from Technical University of Munich, Germany for improving the written quality of this paper.

7. Conclusions

It is shown in Section 4 that to obtain the exact tangent stiffness and tangent matrices of generalized force vectors one must compute such terms ∂ (BV7 )/∂ q and ∂ (BT V3 )/∂ q (superscripts and subscripts for the body index are omitted here for brevity). For the convenience of numerical calculation, in this paragraph concrete expressions of these terms and of the important matrix B are given in a component manner. Define the origin of the local coordinate frame which is measured in the global coordinate frame as R = [x, y, z]T . The orientation of the local coordinate frame described by quaternion is Θ = [e0 ,e1 ,e2 ,e3 ]T . The generalized coordinates of such a strut thus yield q = [RT ,ΘT ]T . Define Vn = [v1 ,v2 ,vn ]T as an arbitrary n dimensional column vector. Define u = [ux ,uy ,uz ]T as the position vector of the associated material point which is measured in the local coordinate frame.

Strut collisions can be frequently encountered in tensegrity static and dynamic analyses and can seriously change the general mechanical behaviors of structures. For many years, there have been few investigations into this topic, which is perhaps due to a lack of effective tools available addressing complex nonlinear effects involved. This paper systematically addresses the issue of strut collisions modeling in tensegrity mechanical analyses. A general and efficient framework for strut collision modeling of tensegrity systems is presented. A two-strut collision model is formulated using contact force methods. Both the normal contact force and the tangential friction force are taken into account. Exact formulations of the generalized force vectors and the related tangent stiffness and tangent damping matrices are analytically derived. Fundamental issues related to collision detection and stiffness coefficient determination are also addressed. Two representative numerical examples are presented to demonstrate the effectiveness of the general framework in addressing strut collision effects. From the first example involving a classical tensegrity unit, we show that the axial stiffness can be significantly enhanced due to strut collisions. The radial dimensions of struts have a direct influence on when this phenomenon begins to take effect. The outcomes presented here may guide the further development of extremal materials or metamaterial designs by taking advantage of such collision effects through the deliberate design of geometry parameters. The second example involves folding analyses of a clustered tensegrity tower. The static results show that by allowing for some strut collisions, one may make full use of structural features to realize, for instance, achieve an extremely low degree of volume occupancy. The dynamic results show that strut collision behavior at the final folded state tends to decrease the amplitude of structural vibration. Total collision forces are dominated by normal forces. In contrast to position curves, collisions forces are very sensitive to the actuation speed. The outcomes thus provide a reference solution for verifying under what control speed the quasi-static solution is applicable. The discussions given here provide us with a basic paradigm on finding a limit actuation speed for deploying/folding such structures. This paper represents a first attempt at strut collision simulation in tensegrity analyses. It should be mentioned that the weakness of the present model is that it is not suited to thin strut analyses for which bending effects may be particularly obvious when collision occurs. In such cases, a flexible model must be used. Compared to the flexible model, the proposed rigid body model may overestimate collision effects. However, this may also be viewed as an advantage of the proposed model, as it renders the design process more conservative, further ensuring structural safety. Future works may be conducted by developing appropriately flexible models and drawing a thorough and quantitative comparison between different mechanical models. Besides, experimental verifications are also of great importance. Acknowledgments The authors are grateful for the financial support of the National Key Research and Development Plan (2016YFB0200702); the National Natural Science Foundation of China (11772074, 11761131005, 91748203, 91515103). The authors would like to thank

Appendix A

A.1. Expressions of associated matrix B ∈ R3×7

B ( 1, 1 ) = 1; B ( 2, 2 ) = 1; B ( 3, 3 ) = 1; B ( 1, 4 ) = 2e0 ux − 2e3 uy + 2e2 uz ; B ( 1, 5 ) = 2e1 ux + 2e2 uy + 2e3 uz ; B ( 1, 6 ) = 2e1 uy − 2e2 ux + 2e0 uz ; B ( 1, 7 ) = 2e1 uz − 2e0 uy − 2e3 ux ; B ( 2, 4 ) = 2e3 ux + 2e0 uy − 2e1 uz ; B ( 2, 5 ) = 2e2 ux − 2e1 uy − 2e0 uz ;

(A.1)

B ( 2, 6 ) = 2e1 ux + 2e2 uy + 2e3 uz ; B ( 2, 7 ) = 2e0 ux − 2e3 uy + 2e2 uz ; B ( 3, 4 ) = 2e1 uy − 2e2 ux + 2e0 uz ; B ( 3, 5 ) = 2e3 ux + 2e0 uy − 2e1 uz ; B ( 3, 6 ) = 2e3 uy − 2e0 ux − 2e2 uz ; B ( 3, 7 ) = 2e1 ux + 2e2 uy + 2e3 uz ; A.2. Expressions of associated matrix M1 = ∂ (BV7 )/∂ q ∈ R3×7

M1 ( 1, 4 ) = 2ux v4 − 2uy v7 + 2uz v6 ; M1 ( 1, 5 ) = 2ux v5 + 2uy v6 + 2uz v7 ; M1 ( 1, 6 ) = 2uy v5 − 2ux v6 + 2uz v4 ; M1 ( 1, 7 ) = 2uz v5 − 2uy v4 − 2ux v7 ; M1 ( 2, 4 ) = 2ux v7 + 2uy v4 − 2uz v5 ; M1 ( 2, 5 ) = 2ux v6 − 2uy v5 − 2uz v4 ; M1 ( 2, 6 ) = 2ux v5 + 2uy v6 + 2uz v7 ; M1 ( 2, 7 ) = 2ux v4 − 2uy v7 + 2uz v6 ; M1 ( 3, 4 ) = 2uy v5 − 2ux v6 + 2uz v4 ; M1 ( 3, 5 ) = 2ux v7 + 2uy v4 − 2uz v5 ; M1 ( 3, 6 ) = 2uy v7 − 2ux v4 − 2uz v6 ; M1 ( 3, 7 ) = 2ux v5 + 2uy v6 + 2uz v7 ;

(A.2)

Z. Kan, H. Peng and B. Chen et al. / International Journal of Solids and Structures 167 (2019) 202–219

A.3. Expressions of associated matrix M2 = ∂ (BT V3 )/∂ q ∈ R7×7

M2 ( 4, 4 ) = 2ux v1 + 2uy v2 + 2uz v3 ; M2 ( 4, 5 ) = 2uy v3 − 2uz v2 ; M2 ( 4, 6 ) = 2uz v1 − 2ux v3 ; M2 ( 4, 7 ) = 2ux v2 − 2uy v1 ; M2 ( 5, 4 ) = 2uy v3 − 2uz v2 ; M2 ( 5, 5 ) = 2ux v1 − 2uy v2 − 2uz v3 ; M2 ( 5, 6 ) = 2ux v2 + 2uy v1 ; M2 ( 5, 7 ) = 2ux v3 + 2uz v1 ; M2 ( 6, 4 ) = 2uz v1 − 2ux v3 ;

(A.3)

M2 ( 6, 5 ) = 2ux v2 + 2uy v1 ; M2 ( 6, 6 ) = 2uy v2 − 2ux v1 − 2uz v3 ; M2 ( 6, 7 ) = 2uy v3 + 2uz v2 ; M2 ( 7, 4 ) = 2ux v2 − 2uy v1 ; M2 ( 7, 5 ) = 2ux v3 + 2uz v1 ; M2 ( 7, 6 ) = 2uy v3 + 2uz v2 ; M2 ( 7, 7 ) = 2uz v3 − 2uy v2 − 2ux v1 ; Components of matrices B, M1 and M2 not presented in Eqs. (A.1), (A.2) and (A.3) are equal to zero. References Acary, V., Brogliato, B., 2008. Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics. Springer-Verlag, Berlin Heidelberg. Al Sabouni-Zawadzka, A., Gilewski, W., 2018. Smart metamaterial based on the simplex tensegrity pattern. Materials 11, 673. Alves, J., Peixinho, N., da Silva, M.T., Flores, P., Lankarani, H.M., 2015. A comparative study of the viscoelastic constitutive models for frictionless contact interfaces in solids. Mech. Mach. Theory 85, 172–188. Ambrósio, J.A.C., 2002. Impact of rigid and flexible multibody systems: deformation description and contact models. In: Schiehlen, W., Valášek, M. (Eds.), Virtual Nonlinear Multibody Systems. Springer, Dordrecht, pp. 57–81. Ashwear, N., Eriksson, A., 2015. Influence of temperature on the vibration properties of tensegrity structures. Int. J. Mech. Sci. 99, 237–250. Barnes, M., 1999. Form finding and analysis of tension structures by dynamic relaxation. Int. J. Space Struct. 14, 89–104. Bel Hadj Ali, N., Rhode-Barbarigos, L., Smith, I.F.C., 2011. Analysis of clustered tensegrity structures using a modified dynamic relaxation algorithm. Int. J. Solids Struct. 48, 637–647. Bottasso, C.L., Dopico, D., Trainelli, L., 2007. On the optimal scaling of index three DAEs in multibody dynamics. Multibody Syst. Dyn. 19, 3–20. Calladine, C.R., 1978. Buckminster Fuller’s “Tensegrity” structures and Clerk Maxwell’s rules for the construction of stiff frames. Int. J. Solids Struct. 14, 161–172. Cefalo, M., Mirats-Tur, J.M., 2011. A comprehensive dynamic model for class-1 tensegrity systems based on quaternions. Int. J. Solids Struct. 48, 785–802. Cefalo, M., Mirats Tur, J.M., 2010. Real-time self-collision detection algorithms for tensegrity systems. Int. J. Solids Struct. 47, 1711–1722. Chen, Q.-z., Acary, V., Virlez, G., Brüls, O., 2013. A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints. Int. J. Numer. Methods Eng. 96, 487–511. Faroughi, S., Khodaparast, H.H., Friswell, M.I., 2015. Non-linear dynamic analysis of tensegrity structures using a co-rotational method. Int. J. Nonlinear Mech. 69, 55–65. Feeny, B., Guran, A.s., Hinrichs, N., Popp, K., 1998. A historical review on dry friction and stick-slip phenomena. Appl. Mech. Rev. 51, 321. Flores, P., Ambrósio, J., Claro, J.C.P., Lankarani, H.M., Koshy, C.S., 2006. A study on dynamics of mechanical systems including joints with clearance and lubrication. Mech. Mach. Theory 41, 247–261.

219

Fraternali, F., Carpentieri, G., Amendola, A., 2015. On the mechanical modeling of the extreme softening/stiffening response of axially loaded tensegrity prisms. J. Mech. Phys. Solids 74, 136–157. Hanaor, A., 2012. Debunking “Tensegrity”—a personal perspective international. J. Space Struct. 27, 179–183. Hertz, H.R., 1881. Über die Berührung fester elastischer Körper. J. Reine und Angew. Math 92, 156–171. Johnson, K.L., 1985. Contact Mechanics. Cambridge University Press, Cambridge, UK. Juan, S.H., Mirats Tur, J.M., 2008. A method to generate stable, collision free configurations for tensegrity based robots. In: 2008 IEEE/RSJ International Conference on Intelligent Robots and Systems. Nice, France, pp. 3769–3774. Juan, S.H., Skelton, R.E., Mirats Tur, J.M., 2009. Dynamically stable collision avoidance for tensegrity based robots. ASME/IFToMM International Conference on Reconfigurable Mechanisms and Robots. Kan, Z., Peng, H., Chen, B., Zhong, W., 2018a. Nonlinear dynamic and deployment analysis of clustered tensegrity structures using a positional formulation FEM. Compos. Struct. 187, 241–258. Kan, Z., Peng, H., Chen, B., Zhong, W., 2018b. A sliding cable element of multibody dynamics with application to nonlinear dynamic deployment analysis of clustered tensegrity. Int. J. Solids Struct. 130–131, 61–79. Kebiche, K., Kazi-Aoual, M.N., Motro, R., 1999. Geometrical non-linear analysis of tensegrity systems. Eng. Struct. 21, 864–876. Korkmaz, S., 2011. A review of active structural control: challenges for engineering informatics. Comput. Struct. 89, 2113–2132. Kwak, B.M., 1991. Complementarity problem formulation of three-dimensional frictional contact. J. Appl. Mech. 58, 134–140. Lankarani, H.M., Nikravesh, P.E., 1990. A contact force model with hysteresis damping for impact analysis of multibody systems. J. Mech. Des. 112, 369. Le Saux, C., Cevaer, F., Motro, R., 2004. Contribution to 3D impact problems: collisions between two slender steel bars. Comptes Rendus Mécanique 332, 17–22. Milton, G.W., Cherkaev, A.V., 1995. Which elasticity tensors are realizable? J. Eng. Mater. Technol. 117, 483. Moored, K.W., Bart-Smith, H., 2009. Investigation of clustered actuation in tensegrity structures. Int. J. Solids Struct. 46, 3272–3281. Murakami, H., 2001. Static and dynamic analyses of tensegrity structures. Part 1. Nonlinear equations of motion. Int. J. Solids Struct. 38, 3599–3613. Nagase, K., Skelton, R.E., 2015. Double-helix tensegrity structures. AIAA J. 53, 847–862. Newmark, N.M., 1959. A method of computation for structural dynamics. J. Eng. Mech. Div. 85, 67–94. Oden, J.T., Martins, J.A.C., 1985. Models and computational methods for dynamic friction phenomena. Comput. Method Appl. Mech. Eng. 52, 527–634. Pellegrino, S., 2001. Deployable Structures. Springer-Verlag Wien, New York. Pfeiffer, F., Glocker, C., 1996. Multibody Dynamics with Unilateral Contacts. John Wiley Sons Inc. Porta, J.M., Hernández-Juan, S., 2016. Path planning for active tensegrity structures. Int. J. Solids Struct. 78-79, 47–56. Schneider, P., Eberly, D., 2002. Geometric Tools for Computer Graphics. Morgan Kaufmann Publishers. Shabana, A.A., 2005. Dynamics of Multibody Systems, Third ed. Cambridge University Press. Skelton, R.E., Oliveira, M.C., 2009. Tensegrity Systems. Springer 2009 ed. Skrinjar, L., Slavicˇ , J., Boltežar, M., 2018. A review of continuous contact-force models in multibody dynamics. Int. J. Mech. Sci. 145, 171–187. Snelson, K., 1965. Continuous tension, discontinuous compression structures. United States Patent 3169611. Sultan, C., 2009. Tensegrity: 60 Years of Art, Science, and Engineering, Adv. Appl. Mech. Elsevier Science & Technology, pp. 60–145. Sultan, C., Corless, M., Skelton, R.E., 20 0 0. Tensegrity flight simulator. J. Guid. Control Dynam. 23, 1055–1064. Tian, Q., Flores, P., Lankarani, H.M., 2018. A comprehensive survey of the analytical, numerical and experimental methodologies for dynamics of multibody mechanical systems with clearance or imperfect joints. Mech. Mach. Theory 122, 1–57. Tran, H.C., Lee, J., 2010. Initial self-stress design of tensegrity grid structures. Comput. Struct. 88, 558–566. Xu, X., Luo, Y., 2012. Collision-free shape control of a plane tensegrity structure using an incremental dynamic relaxation method and a trial-and-error process. In: Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 227, pp. 266–272. Xu, X., Sun, F., Luo, Y., Xu, Y., 2014. Collision-free path planning of tensegrity structures. J. Struct. Eng. 140, 04013084. Yang, S., Sultan, C., 2017. A comparative study on the dynamics of tensegrity-membrane systems based on multiple models. Int. J. Solids Struct. 113-114, 47–69. Zhang, L., Gao, Q., Liu, Y., Zhang, H., 2016. An efficient finite element formulation for nonlinear analysis of clustered tensegrity. Eng. Comput. 33, 252–273. Zhang, L., Lu, M.K., Zhang, H.W., Yan, B., 2015. Geometrically nonlinear elasto-plastic analysis of clustered tensegrity based on the co-rotational approach. Int. J. Mech. Sci. 93, 154–165.