- Email: [email protected]
Soft Machines: Challenges to Computational Dynamics
Available online at www.sciencedirect.com
ScienceDirect Procedia IUTAM 20 (2017) 10 – 17
24th International Congress of Theoretical and Applied Mechanics
Soft Machines: Challenges to Computational Dynamics Haiyan Hu*, Qiang Tian, Cheng Liu School of Aerospace Engineering, Beijing Institute of Technology, Beijing, 100081, China
Abstract The paper surveys the recent studies of authors on the dynamic modeling and simulation of soft machines in the frame of multibody system dynamics. The studies focus on the geometric nonlinearity of coupled overall motion and large deformation of a soft body, the physical nonlinearity of a soft body made of hyper-elastic or elastoplastic materials, the frictional impacts and contacts of soft bodies, and the efficient dynamic computation algorithm of soft multibody systems governed by a set of differential-algebraic equations of very high dimensions. The paper presents the validation of proposed approach through three case studies, including the locomotion of a soft quadrupedal robot, the spinning deployment of a solar sail of spacecraft, and the deployment of a mesh reflector of satellite antenna, as well as the corresponding experimental studies. Finally, the paper gives some remarks on future researches. © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license © 2016 The Authors. Published by Elsevier B.V. Peer-review under responsibility of organizing committee of the 24th (http://creativecommons.org/licenses/by-nc-nd/4.0/). International Congress of Theoretical and Applied Peer-review under responsibility of organizing committeeMechanics of the 24th International Congress of Theoretical and Applied Mechanics Keywords: multibody system dynamics; absolute nodal coordinate formulation; frictional contact; soft robot; deployable space structure
1. Introduction The concept of soft machine covers a wide range of advanced mechanical systems, such as a soft robot handling fragile objects, crawling along tunnel corners or swimming with bio-mimetic fins [1-3], a morphing aircraft [4] and a deployable spacecraft [5]. The soft machine is mainly composed of soft bodies made of soft materials, including polyamide, polyimide and silicon polymer for adapting to complex environments and missions, or dielectric elastomer and ionic polymer metal composite for actuators and artificial muscles [3]. The soft bodies undergo not only large deformations, but also overall motions during the machine operations. Furthermore, the soft bodies may have impacts and contacts with environments due to the machine missions. The dynamic models of soft machines are a kind of soft multibody systems, which give rise to numerous challenges to computational dynamics. They include the geometrical nonlinearities of coupled overall motion and
* Corresponding author. Tel.: +86-10-68915536; fax: +86-10-68912457. E-mail address: [email protected]
2210-9838 © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of organizing committee of the 24th International Congress of Theoretical and Applied Mechanics doi:10.1016/j.piutam.2017.03.003
11
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
large deformation, and the physical nonlinearities due to soft materials and frictional impacts and contacts of soft bodies. Hence, the computational dynamics of soft robots, for example, is far behind the experimental studies of their prototypes [1-3]. This situation definitely blocks the further development of soft robots. As for the morphing aircraft and deployable spacecraft, it is essential to propose proper modeling and efficient computing approaches to the dynamics in their design phase so as to reduce the cost and danger of flight tests [4-5]. The objective of this paper is to present how to deal with the above major challenges of soft machines towards computational dynamics. The rest part of the paper is organized as follows. In Section 2, the dynamic modeling approaches are surveyed via the finite elements of Absolute Nodal Coordinate Formulation (ANCF) and nonlinear constitutive laws of soft materials. Then, the dynamic simulation approaches are outlined with an emphasis on the frictional impacts and contacts of soft bodies in Section 3. In Section 4, three case studies are given to validate the above approaches. Finally, some concluding remarks are made in Section 5. 2. Dynamic Modeling of Soft Bodies 2.1. Description of geometric nonlinearities The geometric nonlinearity of a soft body mainly comes from the coupled overall motion and large deformation of the body. The conventional approaches, such as those based on the floating frame of reference, fail to work for a flexible multibody system when the deformation of any body is large. In recent years, the finite elements of ANCF have served as a powerful tool to describe the coupled overall motion and large deformation of flexible bodies, such as a slender beam [6]. The fully parameterized beam element of ANCF is able to describe the shear and cross-section deformations of the beam, but suffers from locking problems. On the contrary, the gradient deficient beam element of ANCF is suitable for describing the slender beam without shear deformation taken into account and gives good accuracy. In soft machines, the major components are soft bodies, such as cables, membranes and pneumatic chambers. Hence, the gradient deficient finite elements of ANCF play an important role in modelling soft machines. From this consideration, a new spatial curved beam element of 2 nodes (12 nodal coordinates) and a new shell element of 4 nodes (24 nodal coordinates) were established in the frame of gradient deficient ANCF so as to mesh soft bodies [7], and the stiffness reduction criteria were integrated with the new shell element to deal with wrinkles and slacks in membranes [8]. To mesh soft bodies of complicated shapes and guarantee the first-order geometric continuity for adjacent finite elements, several new shell elements of ANCF were also proposed via Bézier triangles and non-uniform rational B-splines [9-10]. The soft machine may also include rigid bodies, which undergo overall motions and also exhibit geometric nonlinearities. As shown in [11], the ANCF reference nodes enable one to describe the rigid bodies via rigid constraint conditions. Hence, it is possible to use a unified approach of ANCF to establish the following set of Differential Algebraic Equations (DAEs) for the dynamics of a soft machine composed of soft and rigid bodies
Md 2 q / dt 2 φqT λ f (q) ® ¯φ(q, t ) 0,
g (q, dq / dt , t ),
(1)
where M R nun is the constant mass matrix, q R n is the generalized coordinate vector, f Rn is the elastic force vector of soft bodies, g Rn is the vector of other forces including friction forces and external forces, φ R m and λ R m (m n) are the constraint vector and the Lagrange multiply vector, respectively. 2.2. Description of physical nonlinearities The soft bodies are usually made of soft materials, such as polymers and elastomers. These soft materials exhibit physical nonlinearities, say, the hyper-elasticity with asymmetric tension and compression behaviors. The
12
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
neu-Hookean constitutive law is a simple model for hyper-elastic materials, and has different versions of compression properties. In a recent study [12], to model soft shells and membranes made of polymers, the gradient deficient shell element of ANCF was integrated with the following version of neu-Hookean constitutive law Q
P I {O[det 2 (C) 1] / 2 P}C1 ,
(2)
where Q is the Piola-Kirchhoff stress tensor of the second kind, C is the right Green strain tensor, O and P are Lame’s constants of hyper-elastic material. The numerical examples of the thin shells subject to large deformation and the inflation of air chambers have well validated the efficacy of the hyper-elastic shell element. Furthermore, the bilinear constitutive law of elastoplastic materials was integrated with the gradient deficient beam element of ANCF [13] to mesh slender beams and thin cables. To reduce computation cost, the bending deformation of the beam element was assumed small and elastic, and only the axial deformation of the beam was elastoplastic. This assumption holds true when the beam element is used to mesh thin cables, which can hardly be subject to any bending loads. 3. Dynamic Simulation of Soft Machines 3.1. Integration of dynamic equations The dynamic simulation of a soft machine is to solve Eq. (1) via a numerical integrator, say, the generalized-α algorithm, step by step in time domain. As Eq. (1) is a set of nonlinear DAEs, the generalized-α algorithm should work with the Newton-Raphson iteration to solve a set of nonlinear algebraic equations at each time step. To reduce computation cost, one can solve linearized algebraic equations in parallel on an OpenMP during the Newton-Raphson iteration [14]. However, it is still an essential issue to efficiently compute the elastic force vector f (q) and the contact force vector in g(q,dq / dt, t ) , which will be addressed in the following two subsections. 3.2. Computation of elastic forces To compute the elastic force vector f e of a fully parameterized finite element of ANCF efficiently, as shown in [14], it is beneficial to rearrange the generalized coordinate vector q e of the finite element as a matrix form q e and factorize it from the integration of strain energy over the element volume v0 such that f e and the corresponding Jacobian yield
fe
K1 (qe )qe , (
wf e )ij wqe
(K1 (q e ))ij (K 2 (q e ))ij ,
1 e 0 0 °(K1 (q ))ij { 2 (qdf qdv qdf qdv ) ³v0 Danbc H ta H kn H fb H vc dv0 , j 3(t 1) m, ® °(K (q e ))ij qmt qdv ³ Danbc H ta H kn H fb H vc dv0 , j 3( f 1) d , v0 ¯ 2
(3)
(4)
where qdf is an entry of q e , Danbc is an entry of constant elasticity matrix D, and H ta is an entry of gradient matrix H of the shape functions, independent of time, of the finite element. Hence, the two integrations in Eq.(4) can be computed and stored before the dynamic simulation. As a result, the computation efficiency of the elastic force vector can be increased about 10 times compared with previous approaches [14]. For the gradient deficient finite elements of ANCF, such as a slender beam element without shear deformation taken into account, it is not possible to factorize all the generalized coordinates from the integration of strain
13
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
energy over the element volume as above. In this case, however, the strain energy of the thin beam includes only axial strain energy and bending strain energy. It is still feasible to factorize the generalized coordinates from the integration of axial strain energy over the element volume and use the simplified computation of bending force [7] so as to speed up the computation of elastic force vector. Regarding to the soft machines of very large scale, it is also possible to reduce some internal degrees of freedom of complicated soft bodies via the sub-structuring techniques in the structure dynamics. This reduction process has proved applicable to some soft bodies discussed in Subsection 4.3. 3.3. Computation of frictional impacts and contacts For the computation of impacts and contacts of two bodies, it is essential to detect when and where the bodies get in touch. To save computation time, it is usual to make the global detection of contacts first by searching a tree of bounding boxes, which cover all possibly contacting bodies, and then to start the detailed detection only when two bounding boxes occupy a same part in space. The detailed detection often begins with the selection of a point on one body, then searches for a corresponding point on the other body so that their distance is minimal, and finally determines the contact when the above distance becomes zero or negative. To detect the contact of two beams of circular section in Fig. 1, it is efficient to skip the global search and to make the detailed search directly since the search along beams is very fast. The two beams are identified as master beam with smaller deformation and slave beam with larger deformation first. Base on selected point N on the central line of slave beam with Rs in radius, the shortest distance from point Q to the master beam with Rm in radius reads [15-16]
g (] s )
min r(] m ) r(] s ) Rm Rs ,
(5)
m
where the pair of points m and s on the central lines of two beams yields the following set of nonlinear algebraic equations T °(r (] m ) r (] s )) r c(] m ) 0, ® T °¯(r (] m ) r (] s )) r c(] s ) 0.
(6)
As shown in Fig. 2, g (] s ) is a continuous function with multiple minima and maxima [16]. The condition g (] s ) d 0 holds true if and only if the two beams get contacts. node midpoint Master beam ric n ri M rc(] m ) t r (] m ) P gm gj rP Q Z rQ o rc(] s ) r (] s ) N X Y rj Slave beam Fig. 1. Distance between two beams
g (] s )
Detection point Extremal point
Extremal value Contact zone
Beam I
Z
g j 1
g j 1
]s
gj g j min( g j 1, g j 1 ) Fig. 2. Distance function along slave beam
X Y
Fn
A Ft
n′ Q′ P′ j' B
i'
D
n i Q P '] j
-Ft
j
C -Fn
Contact zone Beam II
Fig. 3. Regional contact of two beams
Different from relatively rigid bodies, the impacts or contacts of soft bodies are not point-wise, but regional instead, as shown in Fig. 3. Furthermore, the soft body, like a long cable, may get impacts or contacts with itself during an overall motion. Thus, the detection of impacts or contacts of soft bodies becomes more complicated. To deal with the above problems, a detection strategy was proposed for multi-zone contacts of two beams based
14
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
on the shortest distance function g (] s ) [15-16]. The basic idea to determine a contact zone is to locate a pair of contacting points P and Q in Fig. 3 first and then to check whether the adjacent points Pc and Qc keep a contact such that the contact condition g (] j '] j ) d 0 holds true. It is efficient to determine a contact zone of two beams this way after a few iterations. It is intuitive and efficient to compute the normal contact force by using the penalty method even though a proper selection of penalty parameter may need several trials. The tangential contact force is due to a friction process. For soft machines, the relative motion between contacting soft bodies is often slow. In this case, the Coulomb friction model does not offer accurate frictional force including static friction and Stribeck effect, whereas the following LuGre friction model may work well f V z V dz / dt V v , 0 1 2 t °° t z v d / d t f ( v ) z , ® t t ° 1 °¯ f ( v t ) { V 0 v t f n [ Pd ( P s Pd ) exp( v t vs
(7) D
)]1 ,
where ft is the friction force vector, v t is the relative velocity vector, z yields the Ordinary Differential Equation (ODE) in Eq. (7), and other symbols can be found in [16]. The dynamic simulation of two contacting soft bodies needs to solve the ODE at each Gaussian point in a contacting region at each time step. To reduce computation cost, v t is assumed to be unchanged within a fine time step [k 't , (k 1)'t ] such that the ODE in Eq. (7) can be analytically solved [16]. Hence, the following piecewise analytic expression of the LuGre friction model is available for each step of numerical integration
ft( k 1) V 0 z ( k 1) V 1dz ( k 1) / dt V 2 v t( k 1) , ° ( k 1) z ( k ) exp( f ( k 1) 't ) ( f ( k 1) ) 1 v t( k 1) [1 exp( f ( k 1) 't )], ®z ° ( k 1) / dt ( v t( k 1) f ( k 1) z ( k ) ) exp( f ( k 1) 't ), k 1, 2,3, L ¯dz
(8)
4. Case studies of Soft Machines 4.1. Locomotion of a quadrupedal robot The first case study is to demonstrate the validity of the above approach to a simple soft machine, the quadrupedal robot invented by Shepherda et al. [17]. The robot consists of a number of pneumatic chambers made of elastomeric polymers such that the robot is able to move in two gaits, as shown in Fig. 4, under the controlled air pressure in pneumatic chambers. A simplified model for the robot was established via 24 rectangular shell elements of ANCF as shown in Fig. 5. The locomotion of the robot in two gaits was simulated by applying and removing two concentrate forces alternatively in Fig. 5 so as to model the distributed forces of air pressure in pneumatic chambers. Fig. 4 and Fig. 6 present a good agreement between the experimental results and numerical simulations. Furthermore, the numerical simulations indicate that the robot dynamics is robust with respect to the variation of friction coefficients.
Fig. 4. The locomotion of the quadrupedal robot at 5 moments in a test [17].
15
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
F1
F1 T=0s
F2
F2
T=0.6s
T=1.2s
T=1.8s
T=2.4s
Fig. 5. The simplified model and actuation process for the quadrupedal robot in locomotion.
Fig. 6. The locomotion of the quadrupedal robot at 5 moments in a simulation.
4.2. Spinning deployment of a solar sail of spacecraft This case study is to demonstrate the efficacy of the above approach to a practical soft machine, IKAROS launched by Japan Aerospace Exploration Agency in 2010 for deep space exploration. IKAROS is a solar sail composed of 4 trapezoidal polyimide membranes with a thickness of 7μm, and deployed to the square to a length of 14m via the centrifugal forces of 4 lumped masses at the sail corners during the spinning process [5]. To deeply understand the dynamics of the spinning solar sail, a simplified model was numerically and experimentally studied as shown in Fig. 7, where the 4 trapezoidal membranes with the larger base-length of 1.6m were wrapped around a spinning central hub before the test, and bumped out by a releasing device at the initial moment. At the constant rotation speed of 20 rpm of the central hub, the folded membranes were gradually deployed to a square and finally stabilized by the centrifugal forces of the lumped masses at 4 corners of the solar sail. In the numerical study, each membrane was meshed via 16 triangular and 184 rectangular membrane elements of ANCF [18]. In the experimental study, the entire test rig of spinning solar sail was installed in a vacuum chamber to remove the aerodynamic loads on the spinning solar sail. Fig. 7 shows that the numerical results well agree with the experimental results [19]. In addition, a related study showed the importance of proper description of membrane wrinkles. For example, the cyclic “deployment-shrinkage-deployment” phenomenon was observed in numerical simulations if the wrinkles were not properly modeled [8]. A
D
D
D
D
D
C
B
D
A
A
A
A
B
C
0.0s
B
0.1s
C
C
C B
0.25s
B
C A B
0.35s
0.5s
10s
Fig. 7. Dynamic configurations of a spinning solar at 6 specific moments (top subfigures: computed; bottom subfigures: tested).
16
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
4.3. Deployment of a mesh reflector of satellite antenna In order to show the feasibility of the above approach to more complicated soft machines, the deployment of a mesh reflector of satellite antenna was studied as shown in Fig. 8, where the hoop truss to a diameter of 15.6m was modeled via 720 fully parameterized beam elements of ANCF, the front mesh of 648 cables, the rear mesh of 648 cables and 211 adjustable cables were modeled via 32140 gradient deficient beam elements of ANCF, respectively. As a result, the dynamic model of entire reflector includes 199,602 nodal degrees of freedom [20]. To speed up the dynamic simulation, the mesh reflector was decomposed into several independent subsystems by cutting its joint. Then, the Schur complement method was used to eliminate the internal nodal degrees of freedom of subsystem and the Lagrange multipliers for joint constraint equations associated with the internal variables. By using the multilevel decomposition approach, the dimension of the simultaneous linear equations could be further reduced [20]. Fig. 8 presents the dynamic configurations of the mesh reflector during the deployment process at 4 specific moments, respectively. As shown in the Fig. 8, the numerical simulations and experimental results coincide well.
0s
100s
200s
300s
Fig. 8. Dynamic configurations of a mesh reflector at 4 specific moments (left subfigures: computed, right subfigures: tested).
5. Concluding Remarks Soft machine is a new concept for the mechanical systems mainly composed of soft bodies for special missions, such as handling fragile objects, adapting to complicated environments and deploying to desirable configurations. The finite element approach of ANCF, integrated with incremental expressions of nonlinear constitutive law and contact forces, offers a unified and powerful tool for the dynamic modeling of soft machines. The simplification of elastic force vector and its Jacobian of ANCF, together with an OpenMP, makes the above approach applicable to the dynamic simulation of complicated soft machines. The case studies of a soft robot and two deployable space structures well validate the effectiveness of the approach. Compared with an increasing interest in soft robots and the high requirements for delpoyable space structures, the dynamic modeling and simulation have not yet received enough attention. The future researches in the computational dynamics of soft machines in the frame of finite elements of ANCF should cover, but are not limited to, the accurate description of complex geometry of soft bodies via iso-geometric analysis, the efficient computation of restoring force vector and its Jacobian integrated with the constitutive law of new soft materials including dielectric elastomer and ionic polymer metal composite, the efficient computation of frictional impacts and contacts of three dimensional soft bodies, the dynamic analysis of soft machines with uncertainties in
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
frictional impacts and contacts, the reduction of dynamic models of soft machines. Furthermore, the integrated studies, such as the corelation between computations and experiments, and the dynamics and control, should also receive more attention.
Acknowledgements This work was supported in part by The National Natural Science Foundation of China under Grants No. 11290150 and 11290151.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]
Rus D, Tolley MT. Design, fabrication and control of soft robots. Nature 2015; 521: 467-475. Iida F, Laschi C. Soft robotics: challenges and perspectives. Procedia Computer Science 2011; 7: 99-102. Li TF, Li GR, Liang YM, Cheng TY, Yang XX, Huang ZL. Review of materials and structures in soft robotics, Chinese Journal of Theoretical and Applied Mechanics 2016; 48: 756-766. Ajaj RM, Beaverstock CS, Friswell MI. Morphing aircraft: the need for a new design philosophy. Aerospace Science and Technology 2016; 49: 154-166. Tsuda Y, Mori O, Funase R, Sawada H, Yamamoto T, Saiki T, el a. Achievement of IKAROS - Japanese deep space solar sail demonstration mission. Acta Astronautica 2013; 82: 183-188. Gerstmayr J, Sugiyama H, Mikkola A. Review on the absolute nodal coordinate formulation for large deformation analysis of multibody systems, Journal of Computational and Nonlinear Dynamics 2013; 8: 031016. Liu C, Tian Q, Hu HY. New spatial curved beam and shell elements of gradient deficient absolute nodal coordinate formulation, Nonlinear Dynamics 2012; 70: 1903-1918. Liu C, Tian Q, Yan D, Hu HY. Dynamic analysis of membrane systems undergoing overall motions, large deformations and wrinkles via thin shell elements of ANCF. Computer Methods in Applied Mechanics and Engineering 2013; 258: 81-95. Chang HJ, Liu C, Tian Q, Hu HY, Mikkola A. Three new triangular shell elements of ANCF represented by Bézier triangles. Multibody System Dynamics 2015; 35: 321-351. Luo K, Liu C, Tian Q, Hu HY. An efficient model reduction method for buckling analyses of thin shells based on IGA. Computer Methods in Applied Mechanics and Engineering 2016; 309: 243-268. Shabana AA. ANCF reference node for multibody system analysis. Proceedings of the Institution of Mechanical Engineers Part K Journal of Multi-body Dynamics 2014; 229: 109-112. Luo K, Liu C, Tian Q, Hu HY. Nonlinear static and dynamic analysis of hyper-elastic thin shells via the absolute nodal coordinate formulation. Nonlinear Dynamics 2016; 85: 949-971. Wang QT, Tian Q, Hu HY. Contact dynamics of elasto-plastic thin beams simulated via absolute nodal coordinate formulation. Acta Mechanica Sinica 2016; 32: 525-534. Liu C, Tian Q, Hu HY. Dynamics of a large scale rigid-flexible multibody system with composite laminated plates. Multibody System Dynamics 2011; 26: 283-305. Wang QT, Tian Q, Hu HY. Dynamic simulation of frictional contacts of thin beams during large overall motions via absolute nodal coordinate formulation. Nonlinear Dynamics 2014; 77: 1411-1425. Wang QT, Tian Q, Hu HY. Dynamic simulation of frictional multi-zone contacts of thin beams. Nonlinear Dynamics 2016; 83: 19191937. Shepherda RF, Ilievskia F, Choia W, Morina SA, Stokesa AA, Mazzeoa AD, et al. Multigait soft robot. Proceedings of National Academy of Sciences 2011; 108: 20400-20403. Zhao J, Tian Q, Hu HY. Deployment dynamics of a simplified spinning “IKAROS” solar sail via absolute coordinate based method. Acta Mechanica Sinica 2013; 29: 132-142. Zhou XJ, Zhou CY, Zhang XX, Hu HY. Ground simulation tests of spinning deployment dynamics of a solar sail. Journal of Vibration Engineering 2015; 28: 175-182. Li P, Liu C, Tian Q, Hu HY, Song YP. Dynamics of a deployable mesh reflector of satellite antenna: parallel computation and deployment simulation. Journal of Computational and Nonlinear Dynamics 2016; 11: 061005.
17
ScienceDirect Procedia IUTAM 20 (2017) 10 – 17
24th International Congress of Theoretical and Applied Mechanics
Soft Machines: Challenges to Computational Dynamics Haiyan Hu*, Qiang Tian, Cheng Liu School of Aerospace Engineering, Beijing Institute of Technology, Beijing, 100081, China
Abstract The paper surveys the recent studies of authors on the dynamic modeling and simulation of soft machines in the frame of multibody system dynamics. The studies focus on the geometric nonlinearity of coupled overall motion and large deformation of a soft body, the physical nonlinearity of a soft body made of hyper-elastic or elastoplastic materials, the frictional impacts and contacts of soft bodies, and the efficient dynamic computation algorithm of soft multibody systems governed by a set of differential-algebraic equations of very high dimensions. The paper presents the validation of proposed approach through three case studies, including the locomotion of a soft quadrupedal robot, the spinning deployment of a solar sail of spacecraft, and the deployment of a mesh reflector of satellite antenna, as well as the corresponding experimental studies. Finally, the paper gives some remarks on future researches. © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license © 2016 The Authors. Published by Elsevier B.V. Peer-review under responsibility of organizing committee of the 24th (http://creativecommons.org/licenses/by-nc-nd/4.0/). International Congress of Theoretical and Applied Peer-review under responsibility of organizing committeeMechanics of the 24th International Congress of Theoretical and Applied Mechanics Keywords: multibody system dynamics; absolute nodal coordinate formulation; frictional contact; soft robot; deployable space structure
1. Introduction The concept of soft machine covers a wide range of advanced mechanical systems, such as a soft robot handling fragile objects, crawling along tunnel corners or swimming with bio-mimetic fins [1-3], a morphing aircraft [4] and a deployable spacecraft [5]. The soft machine is mainly composed of soft bodies made of soft materials, including polyamide, polyimide and silicon polymer for adapting to complex environments and missions, or dielectric elastomer and ionic polymer metal composite for actuators and artificial muscles [3]. The soft bodies undergo not only large deformations, but also overall motions during the machine operations. Furthermore, the soft bodies may have impacts and contacts with environments due to the machine missions. The dynamic models of soft machines are a kind of soft multibody systems, which give rise to numerous challenges to computational dynamics. They include the geometrical nonlinearities of coupled overall motion and
* Corresponding author. Tel.: +86-10-68915536; fax: +86-10-68912457. E-mail address: [email protected]
2210-9838 © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of organizing committee of the 24th International Congress of Theoretical and Applied Mechanics doi:10.1016/j.piutam.2017.03.003
11
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
large deformation, and the physical nonlinearities due to soft materials and frictional impacts and contacts of soft bodies. Hence, the computational dynamics of soft robots, for example, is far behind the experimental studies of their prototypes [1-3]. This situation definitely blocks the further development of soft robots. As for the morphing aircraft and deployable spacecraft, it is essential to propose proper modeling and efficient computing approaches to the dynamics in their design phase so as to reduce the cost and danger of flight tests [4-5]. The objective of this paper is to present how to deal with the above major challenges of soft machines towards computational dynamics. The rest part of the paper is organized as follows. In Section 2, the dynamic modeling approaches are surveyed via the finite elements of Absolute Nodal Coordinate Formulation (ANCF) and nonlinear constitutive laws of soft materials. Then, the dynamic simulation approaches are outlined with an emphasis on the frictional impacts and contacts of soft bodies in Section 3. In Section 4, three case studies are given to validate the above approaches. Finally, some concluding remarks are made in Section 5. 2. Dynamic Modeling of Soft Bodies 2.1. Description of geometric nonlinearities The geometric nonlinearity of a soft body mainly comes from the coupled overall motion and large deformation of the body. The conventional approaches, such as those based on the floating frame of reference, fail to work for a flexible multibody system when the deformation of any body is large. In recent years, the finite elements of ANCF have served as a powerful tool to describe the coupled overall motion and large deformation of flexible bodies, such as a slender beam [6]. The fully parameterized beam element of ANCF is able to describe the shear and cross-section deformations of the beam, but suffers from locking problems. On the contrary, the gradient deficient beam element of ANCF is suitable for describing the slender beam without shear deformation taken into account and gives good accuracy. In soft machines, the major components are soft bodies, such as cables, membranes and pneumatic chambers. Hence, the gradient deficient finite elements of ANCF play an important role in modelling soft machines. From this consideration, a new spatial curved beam element of 2 nodes (12 nodal coordinates) and a new shell element of 4 nodes (24 nodal coordinates) were established in the frame of gradient deficient ANCF so as to mesh soft bodies [7], and the stiffness reduction criteria were integrated with the new shell element to deal with wrinkles and slacks in membranes [8]. To mesh soft bodies of complicated shapes and guarantee the first-order geometric continuity for adjacent finite elements, several new shell elements of ANCF were also proposed via Bézier triangles and non-uniform rational B-splines [9-10]. The soft machine may also include rigid bodies, which undergo overall motions and also exhibit geometric nonlinearities. As shown in [11], the ANCF reference nodes enable one to describe the rigid bodies via rigid constraint conditions. Hence, it is possible to use a unified approach of ANCF to establish the following set of Differential Algebraic Equations (DAEs) for the dynamics of a soft machine composed of soft and rigid bodies
Md 2 q / dt 2 φqT λ f (q) ® ¯φ(q, t ) 0,
g (q, dq / dt , t ),
(1)
where M R nun is the constant mass matrix, q R n is the generalized coordinate vector, f Rn is the elastic force vector of soft bodies, g Rn is the vector of other forces including friction forces and external forces, φ R m and λ R m (m n) are the constraint vector and the Lagrange multiply vector, respectively. 2.2. Description of physical nonlinearities The soft bodies are usually made of soft materials, such as polymers and elastomers. These soft materials exhibit physical nonlinearities, say, the hyper-elasticity with asymmetric tension and compression behaviors. The
12
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
neu-Hookean constitutive law is a simple model for hyper-elastic materials, and has different versions of compression properties. In a recent study [12], to model soft shells and membranes made of polymers, the gradient deficient shell element of ANCF was integrated with the following version of neu-Hookean constitutive law Q
P I {O[det 2 (C) 1] / 2 P}C1 ,
(2)
where Q is the Piola-Kirchhoff stress tensor of the second kind, C is the right Green strain tensor, O and P are Lame’s constants of hyper-elastic material. The numerical examples of the thin shells subject to large deformation and the inflation of air chambers have well validated the efficacy of the hyper-elastic shell element. Furthermore, the bilinear constitutive law of elastoplastic materials was integrated with the gradient deficient beam element of ANCF [13] to mesh slender beams and thin cables. To reduce computation cost, the bending deformation of the beam element was assumed small and elastic, and only the axial deformation of the beam was elastoplastic. This assumption holds true when the beam element is used to mesh thin cables, which can hardly be subject to any bending loads. 3. Dynamic Simulation of Soft Machines 3.1. Integration of dynamic equations The dynamic simulation of a soft machine is to solve Eq. (1) via a numerical integrator, say, the generalized-α algorithm, step by step in time domain. As Eq. (1) is a set of nonlinear DAEs, the generalized-α algorithm should work with the Newton-Raphson iteration to solve a set of nonlinear algebraic equations at each time step. To reduce computation cost, one can solve linearized algebraic equations in parallel on an OpenMP during the Newton-Raphson iteration [14]. However, it is still an essential issue to efficiently compute the elastic force vector f (q) and the contact force vector in g(q,dq / dt, t ) , which will be addressed in the following two subsections. 3.2. Computation of elastic forces To compute the elastic force vector f e of a fully parameterized finite element of ANCF efficiently, as shown in [14], it is beneficial to rearrange the generalized coordinate vector q e of the finite element as a matrix form q e and factorize it from the integration of strain energy over the element volume v0 such that f e and the corresponding Jacobian yield
fe
K1 (qe )qe , (
wf e )ij wqe
(K1 (q e ))ij (K 2 (q e ))ij ,
1 e 0 0 °(K1 (q ))ij { 2 (qdf qdv qdf qdv ) ³v0 Danbc H ta H kn H fb H vc dv0 , j 3(t 1) m, ® °(K (q e ))ij qmt qdv ³ Danbc H ta H kn H fb H vc dv0 , j 3( f 1) d , v0 ¯ 2
(3)
(4)
where qdf is an entry of q e , Danbc is an entry of constant elasticity matrix D, and H ta is an entry of gradient matrix H of the shape functions, independent of time, of the finite element. Hence, the two integrations in Eq.(4) can be computed and stored before the dynamic simulation. As a result, the computation efficiency of the elastic force vector can be increased about 10 times compared with previous approaches [14]. For the gradient deficient finite elements of ANCF, such as a slender beam element without shear deformation taken into account, it is not possible to factorize all the generalized coordinates from the integration of strain
13
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
energy over the element volume as above. In this case, however, the strain energy of the thin beam includes only axial strain energy and bending strain energy. It is still feasible to factorize the generalized coordinates from the integration of axial strain energy over the element volume and use the simplified computation of bending force [7] so as to speed up the computation of elastic force vector. Regarding to the soft machines of very large scale, it is also possible to reduce some internal degrees of freedom of complicated soft bodies via the sub-structuring techniques in the structure dynamics. This reduction process has proved applicable to some soft bodies discussed in Subsection 4.3. 3.3. Computation of frictional impacts and contacts For the computation of impacts and contacts of two bodies, it is essential to detect when and where the bodies get in touch. To save computation time, it is usual to make the global detection of contacts first by searching a tree of bounding boxes, which cover all possibly contacting bodies, and then to start the detailed detection only when two bounding boxes occupy a same part in space. The detailed detection often begins with the selection of a point on one body, then searches for a corresponding point on the other body so that their distance is minimal, and finally determines the contact when the above distance becomes zero or negative. To detect the contact of two beams of circular section in Fig. 1, it is efficient to skip the global search and to make the detailed search directly since the search along beams is very fast. The two beams are identified as master beam with smaller deformation and slave beam with larger deformation first. Base on selected point N on the central line of slave beam with Rs in radius, the shortest distance from point Q to the master beam with Rm in radius reads [15-16]
g (] s )
min r(] m ) r(] s ) Rm Rs ,
(5)
m
where the pair of points m and s on the central lines of two beams yields the following set of nonlinear algebraic equations T °(r (] m ) r (] s )) r c(] m ) 0, ® T °¯(r (] m ) r (] s )) r c(] s ) 0.
(6)
As shown in Fig. 2, g (] s ) is a continuous function with multiple minima and maxima [16]. The condition g (] s ) d 0 holds true if and only if the two beams get contacts. node midpoint Master beam ric n ri M rc(] m ) t r (] m ) P gm gj rP Q Z rQ o rc(] s ) r (] s ) N X Y rj Slave beam Fig. 1. Distance between two beams
g (] s )
Detection point Extremal point
Extremal value Contact zone
Beam I
Z
g j 1
g j 1
]s
gj g j min( g j 1, g j 1 ) Fig. 2. Distance function along slave beam
X Y
Fn
A Ft
n′ Q′ P′ j' B
i'
D
n i Q P '] j
-Ft
j
C -Fn
Contact zone Beam II
Fig. 3. Regional contact of two beams
Different from relatively rigid bodies, the impacts or contacts of soft bodies are not point-wise, but regional instead, as shown in Fig. 3. Furthermore, the soft body, like a long cable, may get impacts or contacts with itself during an overall motion. Thus, the detection of impacts or contacts of soft bodies becomes more complicated. To deal with the above problems, a detection strategy was proposed for multi-zone contacts of two beams based
14
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
on the shortest distance function g (] s ) [15-16]. The basic idea to determine a contact zone is to locate a pair of contacting points P and Q in Fig. 3 first and then to check whether the adjacent points Pc and Qc keep a contact such that the contact condition g (] j '] j ) d 0 holds true. It is efficient to determine a contact zone of two beams this way after a few iterations. It is intuitive and efficient to compute the normal contact force by using the penalty method even though a proper selection of penalty parameter may need several trials. The tangential contact force is due to a friction process. For soft machines, the relative motion between contacting soft bodies is often slow. In this case, the Coulomb friction model does not offer accurate frictional force including static friction and Stribeck effect, whereas the following LuGre friction model may work well f V z V dz / dt V v , 0 1 2 t °° t z v d / d t f ( v ) z , ® t t ° 1 °¯ f ( v t ) { V 0 v t f n [ Pd ( P s Pd ) exp( v t vs
(7) D
)]1 ,
where ft is the friction force vector, v t is the relative velocity vector, z yields the Ordinary Differential Equation (ODE) in Eq. (7), and other symbols can be found in [16]. The dynamic simulation of two contacting soft bodies needs to solve the ODE at each Gaussian point in a contacting region at each time step. To reduce computation cost, v t is assumed to be unchanged within a fine time step [k 't , (k 1)'t ] such that the ODE in Eq. (7) can be analytically solved [16]. Hence, the following piecewise analytic expression of the LuGre friction model is available for each step of numerical integration
ft( k 1) V 0 z ( k 1) V 1dz ( k 1) / dt V 2 v t( k 1) , ° ( k 1) z ( k ) exp( f ( k 1) 't ) ( f ( k 1) ) 1 v t( k 1) [1 exp( f ( k 1) 't )], ®z ° ( k 1) / dt ( v t( k 1) f ( k 1) z ( k ) ) exp( f ( k 1) 't ), k 1, 2,3, L ¯dz
(8)
4. Case studies of Soft Machines 4.1. Locomotion of a quadrupedal robot The first case study is to demonstrate the validity of the above approach to a simple soft machine, the quadrupedal robot invented by Shepherda et al. [17]. The robot consists of a number of pneumatic chambers made of elastomeric polymers such that the robot is able to move in two gaits, as shown in Fig. 4, under the controlled air pressure in pneumatic chambers. A simplified model for the robot was established via 24 rectangular shell elements of ANCF as shown in Fig. 5. The locomotion of the robot in two gaits was simulated by applying and removing two concentrate forces alternatively in Fig. 5 so as to model the distributed forces of air pressure in pneumatic chambers. Fig. 4 and Fig. 6 present a good agreement between the experimental results and numerical simulations. Furthermore, the numerical simulations indicate that the robot dynamics is robust with respect to the variation of friction coefficients.
Fig. 4. The locomotion of the quadrupedal robot at 5 moments in a test [17].
15
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
F1
F1 T=0s
F2
F2
T=0.6s
T=1.2s
T=1.8s
T=2.4s
Fig. 5. The simplified model and actuation process for the quadrupedal robot in locomotion.
Fig. 6. The locomotion of the quadrupedal robot at 5 moments in a simulation.
4.2. Spinning deployment of a solar sail of spacecraft This case study is to demonstrate the efficacy of the above approach to a practical soft machine, IKAROS launched by Japan Aerospace Exploration Agency in 2010 for deep space exploration. IKAROS is a solar sail composed of 4 trapezoidal polyimide membranes with a thickness of 7μm, and deployed to the square to a length of 14m via the centrifugal forces of 4 lumped masses at the sail corners during the spinning process [5]. To deeply understand the dynamics of the spinning solar sail, a simplified model was numerically and experimentally studied as shown in Fig. 7, where the 4 trapezoidal membranes with the larger base-length of 1.6m were wrapped around a spinning central hub before the test, and bumped out by a releasing device at the initial moment. At the constant rotation speed of 20 rpm of the central hub, the folded membranes were gradually deployed to a square and finally stabilized by the centrifugal forces of the lumped masses at 4 corners of the solar sail. In the numerical study, each membrane was meshed via 16 triangular and 184 rectangular membrane elements of ANCF [18]. In the experimental study, the entire test rig of spinning solar sail was installed in a vacuum chamber to remove the aerodynamic loads on the spinning solar sail. Fig. 7 shows that the numerical results well agree with the experimental results [19]. In addition, a related study showed the importance of proper description of membrane wrinkles. For example, the cyclic “deployment-shrinkage-deployment” phenomenon was observed in numerical simulations if the wrinkles were not properly modeled [8]. A
D
D
D
D
D
C
B
D
A
A
A
A
B
C
0.0s
B
0.1s
C
C
C B
0.25s
B
C A B
0.35s
0.5s
10s
Fig. 7. Dynamic configurations of a spinning solar at 6 specific moments (top subfigures: computed; bottom subfigures: tested).
16
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
4.3. Deployment of a mesh reflector of satellite antenna In order to show the feasibility of the above approach to more complicated soft machines, the deployment of a mesh reflector of satellite antenna was studied as shown in Fig. 8, where the hoop truss to a diameter of 15.6m was modeled via 720 fully parameterized beam elements of ANCF, the front mesh of 648 cables, the rear mesh of 648 cables and 211 adjustable cables were modeled via 32140 gradient deficient beam elements of ANCF, respectively. As a result, the dynamic model of entire reflector includes 199,602 nodal degrees of freedom [20]. To speed up the dynamic simulation, the mesh reflector was decomposed into several independent subsystems by cutting its joint. Then, the Schur complement method was used to eliminate the internal nodal degrees of freedom of subsystem and the Lagrange multipliers for joint constraint equations associated with the internal variables. By using the multilevel decomposition approach, the dimension of the simultaneous linear equations could be further reduced [20]. Fig. 8 presents the dynamic configurations of the mesh reflector during the deployment process at 4 specific moments, respectively. As shown in the Fig. 8, the numerical simulations and experimental results coincide well.
0s
100s
200s
300s
Fig. 8. Dynamic configurations of a mesh reflector at 4 specific moments (left subfigures: computed, right subfigures: tested).
5. Concluding Remarks Soft machine is a new concept for the mechanical systems mainly composed of soft bodies for special missions, such as handling fragile objects, adapting to complicated environments and deploying to desirable configurations. The finite element approach of ANCF, integrated with incremental expressions of nonlinear constitutive law and contact forces, offers a unified and powerful tool for the dynamic modeling of soft machines. The simplification of elastic force vector and its Jacobian of ANCF, together with an OpenMP, makes the above approach applicable to the dynamic simulation of complicated soft machines. The case studies of a soft robot and two deployable space structures well validate the effectiveness of the approach. Compared with an increasing interest in soft robots and the high requirements for delpoyable space structures, the dynamic modeling and simulation have not yet received enough attention. The future researches in the computational dynamics of soft machines in the frame of finite elements of ANCF should cover, but are not limited to, the accurate description of complex geometry of soft bodies via iso-geometric analysis, the efficient computation of restoring force vector and its Jacobian integrated with the constitutive law of new soft materials including dielectric elastomer and ionic polymer metal composite, the efficient computation of frictional impacts and contacts of three dimensional soft bodies, the dynamic analysis of soft machines with uncertainties in
Haiyan Hu et al. / Procedia IUTAM 20 (2017) 10 – 17
frictional impacts and contacts, the reduction of dynamic models of soft machines. Furthermore, the integrated studies, such as the corelation between computations and experiments, and the dynamics and control, should also receive more attention.
Acknowledgements This work was supported in part by The National Natural Science Foundation of China under Grants No. 11290150 and 11290151.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]
Rus D, Tolley MT. Design, fabrication and control of soft robots. Nature 2015; 521: 467-475. Iida F, Laschi C. Soft robotics: challenges and perspectives. Procedia Computer Science 2011; 7: 99-102. Li TF, Li GR, Liang YM, Cheng TY, Yang XX, Huang ZL. Review of materials and structures in soft robotics, Chinese Journal of Theoretical and Applied Mechanics 2016; 48: 756-766. Ajaj RM, Beaverstock CS, Friswell MI. Morphing aircraft: the need for a new design philosophy. Aerospace Science and Technology 2016; 49: 154-166. Tsuda Y, Mori O, Funase R, Sawada H, Yamamoto T, Saiki T, el a. Achievement of IKAROS - Japanese deep space solar sail demonstration mission. Acta Astronautica 2013; 82: 183-188. Gerstmayr J, Sugiyama H, Mikkola A. Review on the absolute nodal coordinate formulation for large deformation analysis of multibody systems, Journal of Computational and Nonlinear Dynamics 2013; 8: 031016. Liu C, Tian Q, Hu HY. New spatial curved beam and shell elements of gradient deficient absolute nodal coordinate formulation, Nonlinear Dynamics 2012; 70: 1903-1918. Liu C, Tian Q, Yan D, Hu HY. Dynamic analysis of membrane systems undergoing overall motions, large deformations and wrinkles via thin shell elements of ANCF. Computer Methods in Applied Mechanics and Engineering 2013; 258: 81-95. Chang HJ, Liu C, Tian Q, Hu HY, Mikkola A. Three new triangular shell elements of ANCF represented by Bézier triangles. Multibody System Dynamics 2015; 35: 321-351. Luo K, Liu C, Tian Q, Hu HY. An efficient model reduction method for buckling analyses of thin shells based on IGA. Computer Methods in Applied Mechanics and Engineering 2016; 309: 243-268. Shabana AA. ANCF reference node for multibody system analysis. Proceedings of the Institution of Mechanical Engineers Part K Journal of Multi-body Dynamics 2014; 229: 109-112. Luo K, Liu C, Tian Q, Hu HY. Nonlinear static and dynamic analysis of hyper-elastic thin shells via the absolute nodal coordinate formulation. Nonlinear Dynamics 2016; 85: 949-971. Wang QT, Tian Q, Hu HY. Contact dynamics of elasto-plastic thin beams simulated via absolute nodal coordinate formulation. Acta Mechanica Sinica 2016; 32: 525-534. Liu C, Tian Q, Hu HY. Dynamics of a large scale rigid-flexible multibody system with composite laminated plates. Multibody System Dynamics 2011; 26: 283-305. Wang QT, Tian Q, Hu HY. Dynamic simulation of frictional contacts of thin beams during large overall motions via absolute nodal coordinate formulation. Nonlinear Dynamics 2014; 77: 1411-1425. Wang QT, Tian Q, Hu HY. Dynamic simulation of frictional multi-zone contacts of thin beams. Nonlinear Dynamics 2016; 83: 19191937. Shepherda RF, Ilievskia F, Choia W, Morina SA, Stokesa AA, Mazzeoa AD, et al. Multigait soft robot. Proceedings of National Academy of Sciences 2011; 108: 20400-20403. Zhao J, Tian Q, Hu HY. Deployment dynamics of a simplified spinning “IKAROS” solar sail via absolute coordinate based method. Acta Mechanica Sinica 2013; 29: 132-142. Zhou XJ, Zhou CY, Zhang XX, Hu HY. Ground simulation tests of spinning deployment dynamics of a solar sail. Journal of Vibration Engineering 2015; 28: 175-182. Li P, Liu C, Tian Q, Hu HY, Song YP. Dynamics of a deployable mesh reflector of satellite antenna: parallel computation and deployment simulation. Journal of Computational and Nonlinear Dynamics 2016; 11: 061005.
17