Analysis of contact-impact dynamics of soft finger tapping system by using hybrid computational model

Accepted Manuscript

Analysis of contact–impact dynamics of soft finger tapping system by using hybrid computational model Jiongcan Yang , Yunian Shen PII: DOI: Reference:

S0307-904X(19)30210-0 https://doi.org/10.1016/j.apm.2019.04.020 APM 12767

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

30 October 2018 29 March 2019 4 April 2019

Please cite this article as: Jiongcan Yang , Yunian Shen , Analysis of contact–impact dynamics of soft finger tapping system by using hybrid computational model, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.04.020

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights 

An efficient hybrid computational model to analyze the frictional impact responses during tapping of soft finger is developed.



A lumped–parameter model taking the tangential compliance into account is presented to calculate contact forces. Event driven scheme and computational strategy are given to solve non-smooth

CR IP T



differential equation. 

Numerical results are compared with three-dimensional LS-DYNA solutions.



A ‘reserve slip’ phenomenon caused by large tangential contact compliance is

AC

CE

PT

ED

M

AN US

captured.

1

ACCEPTED MANUSCRIPT

Analysis of contact–impact dynamics of soft finger tapping system by using hybrid computational model Jiongcan Yang, Yunian Shen Department of Mechanics and Engineering Science, School of science, Nanjing University of Science and Technology, Nanjing, 210094, P. R. China

CR IP T

Abstract

CE

PT

ED

M

AN US

The aim of this paper is to develop an efficient hybrid computational model to analyze the frictional impact dynamic responses of soft finger during the tapping event. The large deformation field and inertial field of soft structure are discretized by absolute nodal coordinate formulation. Lagrange multiplier method is adopted to account for the constraints between neighbor phalanges. Considering tangential contact compliance, a lumped–parameter model at the local contact zone is presented to calculate contact forces. The governing equations of soft finger tapping system expressed by generalized coordinates are derived. An event driven scheme is given to analyze and evaluate the stick–slip transitions. The governing equation is integrated by generalized–α method. The applications of the hybrid computational model are demonstrated using various soft finger tapping systems acted by different driven moments. The feasibility of the proposed model is validated by comparing with the LS–DYNA solution. The error of the solution calculated by the proposed model for the peak value of contact force is less than 8%. Furthermore, numerical results show that the large structural compliance and driven moments have significant effect on the frictional impact responses. The normal relative motion between the fingertip and the rough target surface will experience 1~4 compression–restitution transitions when Young’s Modulus is from 0.01GPa to 1GPa or the slenderness ratio of phalanx is from 10.7 to 32. When the posture of soft finger is convex at impact instance, the tangential relative motion will experience 3 slip–stick transitions during the contact process. In addition, it also can be found that the tangential contact compliance can reverse the direction of slip (i.e. ‘reverse slip’ phenomenon).

Key words: hybrid computational model; finger; frictional impact; absolute nodal

AC

coordinate formulation; stick–slip; soft robot.

Corresponding author. Tel.: +86-25-84315590; Fax: +86-25-84315875. E-mail addresses: [email protected] or [email protected] (Yunian Shen). 2

ACCEPTED MANUSCRIPT

1. Introduction

AC

CE

PT

ED

M

AN US

CR IP T

Soft robotic finger, which is made almost entirely by soft material, is suitable for the applications in uncertain and dynamic task environments [1]. Due to the intrinsic softness and flexibility, soft finger has some advantages such as great versatility, strong shock absorption ability and good human–computer interaction [2]. Soft robots are beginning to spring up recently, and they have been preliminarily designed in various engineering fields including medicine [3], biology [4], machinery [5] and etc. The frictional impact will happen inevitably when soft finger conducts a mission [6]. For a tapping event of soft finger, it is generally accompanied with some typical phenomena of frictional impact such as impact vibration, structural large deformation, transient wave propagation and stress concentration at the local contact zone [2, 7]. These characteristics will lead to many adverse effects for the system, e.g. the decrease of working accuracy, the increase of noise, even the failure of structure. For the frictional impact problem of soft material, the structure that is in contact does not only have a unilateral constraint in the normal direction of contact surface, but also a frictional constraint in the tangential direction [8–10]. The addition and deletion of the contact constraint will greatly change the vibration mode of the deformable body [11]. Comparing with collinear impact, the transient characteristics of frictional impact are more complicated. That is because there is stick–slip motion (sometimes the slip will reverse) at the contact point [10]. Johnson [12] and Stronge [13] investigate a singular phenomenon that the horizontal translation velocity and the angular velocity of soft ‘super ball’ will reverse direction (see Section 5.2.6 in reference [13]) after it oblique impacts against a rough surface. They conclude that the large tangential contact compliance of soft material is the root reason of this phenomenon and it cannot be explained by rigid impact model. Pfeiffer [14] experimentally investigates the energy variation of the system for a rubber disc obliquely impacting against a half–space. The tangential velocity reversibility is also obtained. To describe and simulate this singular phenomenon very well. Shen and Stronge [15] propose a lumped–parameter contact model, which can analyze the reverse sliding, to discuss the effect of the tangential contact compliance on frictional responses including the contact forces and the tangential relative velocities between neighbor contact points. However, in the current study, only the compliance of local contact zone is taken into account, the rest of the structure is still treated as rigid body. Besides the treatment of contact constraint at the local contact zone, it also needs a suitable dynamic modelling method to calculate the transient responses of the structure locating at the non–contact zone. The early soft finger is usually designed as a multi–link hard body with a soft fingertip which will contact with target body [16]. This type of soft finger is called as ‘soft tip finger’. Generally, for soft tip finger, the phalanx is assumed as rigid body and the transient responses of the system is obtained based on rigid body dynamic theory combined with a local compliant contact model [17–18]. For example, Doulgeri and Arimoto [19] simplify the soft tip finger into a three–degree–of–freedom rigid link with contact damping forces and torques, then the 3

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

large overall motion and impact responses are calculated by conventional rigid body dynamic theory. Recently, the phalanges of soft finger are entirely made by soft material, the dynamic modelling of soft finger faces the difficulties of describing the large dynamic deformation and large overall motion of soft structure at the same time. The rigid body dynamic theory is no longer applicable to the novel soft finger. For resolving this problem, most scholars generally recourse to commercial finite element method (FEM). For example, Yuk and Lin [20] use ABAQUS software to simulate the dynamic behavior of a single–phalanx soft structure made by hydrogel material. The large bending deformation and overall motion is studied when it kicks a rubber ball or grasping a fish. Chen and Zhang [21] also use ABAQUS software to calculate the stress and strain distribution of a single–phalanx soft finger made by thermoplastic elastomer (TPE), which is under normal contact force and tangential frictional force. Devi and Udupa [22] develop a multi–body soft finger, which is more similar to the human beings' actual cases. It is composed of two hard phalanges and two soft phalanges. Its contact event is also analyzed by commercial FEM. However, the conventional commercial FEM is time–consuming [23] because generally a dense refined spatial discretization and a small size of time step [24–25] are necessary to satisfy the contact algorithm, such as penalty function method or Lagrange multiplier method. Obviously, if the soft finger has a large overall motion for a long time before the contact–impact, the conventional commercial FEM is inefficient to calculate the displacement and velocity. Hence, to quickly and systematically analyze the influence of system parameters on transient responses, it is urgent to develop a new and efficient modelling approach. Cha and Koh [16] establish a multi–body model of a humanoid soft finger. Then, the static deflection of the soft finger is obtained through the equilibrium equations and the relationship of deformation geometry. To analyze dynamic problem, the inertial must be considered. Therefore, an appropriate approach should be presented to discretize the structural deformation field and inertial field. The absolute nodal coordinate formulation (ANCF) is an optional method, which can describe the large dynamic deformation well. Hu and Tian [26] successfully use it to discretize the soft structure when studying the contact walking dynamics of soft quadrupedal robot. Their results show ANCF might have great potential to be applied to analyze the contact–impact event of soft finger. In this paper, considering both tangential and normal contact compliance, a hybrid computational model (HCM) is developed to analyze the frictional impact for a humanoid soft finger tapping against a rough target surface (see Fig. 1). Based on the proposed HCM, the structural stress, large deflection and contact forces during the contact–impact phase are calculated. And the normal and tangential relative motions (i.e. compression–restitution and stick–slip) between the contact points are depicted. In Section 2, the HCM is developed. The humanoid soft finger consisting of multiple phalanges and joints is modelled as a flexible multi–link in the HCM. The model is not hard to be degraded for solving the most developed soft fingers consisting of single phalanx with joint. The ANCF is used to discretize the deformation field and inertia filed of soft structure. A lumped–parameter model (LPM) is presented to 4

ACCEPTED MANUSCRIPT

CR IP T

account for the contact constraints. Be in comparison with the penalty function method adopted in LS–DYNA software, the LPM needn’t numerical iteration which will cost tremendous computational time. The governing equations of the whole system expressed by generalized coordinates are derived. In Section 3, the event driven scheme and the computational strategy to analyze stick–slip transitions are given. In Section 4, the comparison between the LS–DYNA solution and the HCM solution is used to validate the feasibility of the HCM. The influence of coefficient of friction, structural compliance and driven moments on the transient responses is discussed.

AN US

Support

Rough target surface

phalanx fixed hinge movable hinge

Fig. 1. Soft robotic finger taps rough target surface

M

2. Hybrid computational model and governing equations

AC

CE

PT

ED

Here, a HCM is developed to study the frictional impact of the soft finger tapping system (see Fig. 2). The HCM is composed of a flexible multi–link combined with a LPM. The flexible multi–link consists of three links with a half–spherical contact end. Each link w has Young’s Modulus Ew, density ρw, diameter dw, Poisson’s ratio σw and length Lw. There is a fixed hinge at the point O, and the neighbour phalanges are connected by a movable hinge (see Fig. 2a). The driven moments are acted on the hinges. To account for the contact constraints, a LPM at the local contact zone is proposed (see Fig. 2b). The local deformation is represented by compliant elements in both tangential and normal direction. These compliant elements are connected to a massless particle P which can either stick or slide on the rough target surface. The HCM has the following assumptions: 1) All components make the motion in the X– O–Y plane and there is no rotation around the axial of the links; 2) The links are seen as Euler–Bernoulli beams, i.e. the effect of shear deformation is neglected; 3) The cross–sections of the beams keep as plane and perpendicular to the axial; 4) The material is isotropic and linear elastic as well as no plastic deformation occurs during the contact–impact.

5

ACCEPTED MANUSCRIPT

Y initial position (a)

τ3

o

τ2

Y

X

(b)

(c)

o

τ2 H rough target Fn C surface

τ1

τ1

Ft

τ3

P Fn

τ2

v C

Cʹ

τ3

node j node k

rough target surface

Ft

node 1 X

Cʹ

fixed hinge movable hinge

Fig. 2. Hybrid computational model for soft finger tapping system with frictional impact: (a)

CR IP T

Flexible multi–link; (b) LPM; (c) Discretization of flexible multi–link.

AC

CE

PT

ED

M

AN US

Considering the transitions of contact topology, a single tapping event should be assumed to consist of three phases: (i) Pre–impact. The soft finger is in large overall spatial motion acted by gravity and driven moments. There is no active constraint between the fingertip and the rough target surface in this phase. (ii) Contact–impact. Frictional impact will occur once the gap between the fingertip and the rough target surface is equal to zero. There is deformation in both normal and tangential direction of the local contact zone. There are a pair of forces, whose value are equivalent and direction are opposite, acting on the fingertip and the rough target surface, respectively. The impact loading will lead to the axial stretch–compression deformations as well as the bending deformations. (iii) Post–impact. The fingertip is out of contact. It will fly away from the rough target surface. But due to the gravity or driven moments, the fingertip may approach the rough target surface again in some cases. Because of the addition–deletion of contact constraint, different phases have respective governing equations and initial conditions. The displacement and velocity distributions in the phalanges at the termination of the contact–impact phase are used as the initial value conditions of the governing equation during the post–impact phase. The geometric nonlinearity of the soft finger mainly comes from the coupling of overall motion and large structural deformation. The conventional approaches, such as those based on the floating frame of reference, are failed to analyze the soft finger tapping event when the deformation of phalanx is large. At present, the finite element method based on ANCF have served as a powerful tool to describe the coupling effect of overall motion and large structural deformation. The fully parameterized beam element of ANCF is able to describe the shear deformations of the beam, but often suffers from locking problems. On the contrary, the gradient deficient beam element of ANCF is suitable for describing the bending deformation of the slender beam without shear deformation and gives good accuracy [25]. Here, we adopt the gradient deficient beam element (i.e. one–dimensional two–node beam element) to discretize the finger structure (see Fig. 2a). In the discretization model (see Fig. 2c), the wth (w=1, 2 and 3) phalanx is divided into nw elements. The total element number of soft finger is

6

ACCEPTED MANUSCRIPT

3

N   nw w1

During the contact–impact phase, the contact forces are calculated by the LPM, which is acted on the node at the fingertip. 2.1. Governing equation of tapping system during non–impact phase

x

node k

undeformed element node j

deformed element q7 k ζ q8 k

AN US

o q4

M

Y

j

ζj q3

rk

q6

CR IP T

For the pre–impact and post–impact phases, the soft finger is in large overall spatial motion acted by gravity and driven moments. There is no contact constraint between the fingertip and the rough target surface. The governing equations of the soft finger have the same expression for both non–impact phases.

q2

q5

rj

X

q1 O

ED

Fig. 3. Global positions of nodes on the sth beam element.

2.1.1. Governing equations of single beam element

CE

PT

The global positions of nodes on the sth (s=1, 2…, N) beam element are illustrated in Fig. 3. The displacement field of a one–dimensional beam element can be obtained using the following polynomial expression:

r

(s)

 rX( s )   a0  a1 x  a2 x 2  a3 x3    (s)    2 3   rY   b0  b1 x  b2 x  b3 x 

(1)

AC

where r ( s ) is the global position vector of an arbitrary point on the sth element. rX( s ) and rY( s ) are the components of r ( s ) in the X and Y direction, respectively. x is the local coordinate of an arbitrary point along the axis of the sth element in the undeformed configuration (0≤x≤l). l is the original length of the beam element (at node j, x=0; while at node k, x=l). Furthermore, the global position vector r ( s ) is defined in terms of the nodal coordinates qe( s ) and the element shape function S. It can be written in a matrix form as follows [27]: r ( s )  Sqe( s )

(2)

The expression of qe( s ) is 7

ACCEPTED MANUSCRIPT

qe( s )   q1 q2

q3

q4

q5

q6

q8 

T

q7

This absolute nodal coordinates in qe( s ) includes the global displacements of the element nodes, which are

q1  rX( s )

x 0

, q2  rY( s )

x 0

, q5  rX( s )

x l

, q6  rY( s )

x l

and the global slopes of the element nodes, which are

rX( s ) x

, q4  x 0

rY( s ) x

, q7  x 0

rX( s ) x

, q8  x l

rY( s ) x

x l

CR IP T

q3 

According to Eq. 2, the global shape function S can be written as

s S 1 0

0 s1

where the functions si=si(ς) ( i  1,

s2 0

0 s2

s3 0

0 s3

s4 0

, 4 ) are defined as

0 s4 

(3)

AN US

s1  1  3 2  2 3 , s2  l   2 2   3  , s3  3 2  2 3 , s4  l  3   2  and ς=x/l. The kinetic energy Te( s ) of the element s can be expressed as 1 1 Te( s )    ( s ) r ( s )T r ( s ) dV  qe( s )T M e( s ) qe( s ) V 2 2 (s) where M e is the constant mass matrix, which can be written as

(4)

l

M e( s )    ( s ) A( s ) S T Sdx

M

0

where  is the density of the element s. The total strain energy U e( s ) of the element s is given as

ED

(s)

1 l  (s) (s) (s) 2 (s) (s) (s) 2  E A   E I       dx a  2 0  1  qe( s )T  K a( s )  K b( s )  qe( s ) 2

PT

U e( s )  U a( s )  U b( s ) 

(5)

AC

CE

where E ( s ) , A( s ) and I ( s ) are Young’s Modulus, cross–sectional area and moment of inertia of the element s, respectively. U a( s ) and U b( s ) are the strain energy caused by axial deformation and bending deformation, respectively. The axial strain  a( s ) and curvature  ( s ) are expressed as

  a( s )   1  q ( s )T S q ( s )  1  e l e  (s) 2    2      qe( s )T S T Sqe( s ) 

where Sl  S T S , S and S are the first and second derivatives of the shape function S with respect to x. Furthermore, the axial stiffness matrix K a( s ) and bending stiffness matrix K b( s ) are as follows, l l l 1 K a( s )  E ( s ) A( s )  Sl qe( s ) qe( s ) Sl dx   Sl dx , K b( s )   E ( s ) I ( s ) S T Sdx 0 0 0 2 (s) The total stiffness matrix K e of the element s is the sum of K a( s ) and K b( s ) , which





8

ACCEPTED MANUSCRIPT

is

K e( s )  K a( s )  K b( s )

(6)

When the soft finger is only acted by gravity force and there is no driven moments at the hinges, according to the principle of virtual work, the generalized external force matrix Qe( s ) of the element s can be written as l 0 Qe( s )    ( s ) A( s ) S T   dx 0 g

(7)

Me( s ) qe( s )  K e( s ) qe( s )  Qe( s )

where qe( s ) is generalized acceleration.

(8)

AN US

2.1.2. Governing equation of single phalanx

CR IP T

where g is the acceleration of gravity. Substituting Eqs. 4, 5 and 7 into the second Lagrange equation, the dynamic equation of the element s can be obtained as

Firstly, assemble all elements of one phalanx without considering the interfacial connection conditions between the neighbor elements, the governing equation of the wth (w=1, 2 and 3) phalanx can be expressed as

where

0 M e(2) 0 0

ED

 M e(1)  0   0   0

M

M p( w) qp( w)  K p( w) qp( w)  Qp( w)

PT

M p( w)

0 0

0

(9)

 qe(1)  0    (2)  0  q ( w) , qp   e    0    ( nw )  M e( nw )  qe 

0 K e(2) 0 0

0 0

0 0 0

AC

CE

  Qe(1)    (2)  Q ( w) ( w)  Kp , Qp   e       (n )  0 K e( nw )  Qe w  In Eq. 9, there are many non–independent or repeated interfacial DOFs. Hence for coupling all elements, a Boolean matrix B obtained by interfacial connection conditions is used to eliminate the repeated DOFs in Eq. 9. Sometimes this process is called the second coordinate transformation, which can be expressed as  K e(1)  0   0   0

qp( w)  Bqp( w) where qp( w) is the generalized displacement vector of the wth phalanx without repeated DOFs. The dynamic equation of a single phalanx after coupling is

M p( w) qp( w)  K p( w) qp( w)  Qp( w)

(10)

where qp( w) , M p( w) , K p( w) and Qp( w) are the generalized acceleration vector, mass 9

ACCEPTED MANUSCRIPT matrix, stiffness matrix and external force vector of the wth phalanx, respectively. The specific expression of the matrixes is

qp( w)  BT qp( w) , M p( w)  BT M p( w) B , K p( w)  BT K p( w) B , Qp( w)  BTQp( w) 2.1.3. Governing equation of whole soft finger By using Lagrange multiplier method to describe the hinge constraint between two neighbor phalanges, the dynamic equation of whole soft finger is written as

Mq  ΦqT λ  Kq  Q

CR IP T

(11)

where Φ is the constraint matrix of whole soft finger. λ is Lagrange multiplier. q and q are the generalized displacement and acceleration vector of whole soft finger, respectively. M, K and Q are the mass matrix, stiffness matrix and generalized external force matrix of whole soft finger, respectively. Their specific expressions are  M p(1)  qp(1)     q  qp(2)  , M   0  0  qp(3)    

 K p(1) 0    0 , K   0  0 M p(3)  

0 K p(2) 0

 Qp(1)  0     0  , Q  Qp(2)   Qp(3)  K p(3)   

AN US

0 M p(2) 0

If the soft finger has driven moments acted on the hinges, the dynamic equation of whole soft finger should include generalized moment matrix τ. Take the sth element as an example, the rotational angle ζj and the generalized coordinates of node j have the following relation:

 sin  j  1  q3  cos  j  f j  q4

q4  q3 

(12)

ED

M

cos  j  sin  j 

PT

where f j  (q32  q42 )1/2 . It is assumed that there exists a driven moment  acting on node j. By the use of the principle of virtual work, the generalized moment matrix τ ( s ) of the element s can be written as

τ ( s )  0 0  q4 / f j2  q3 / f j2

0 0 0 0

T

(13)

CE

In this paper, the driven moments may act on three hinges. Therefore, the generalized moment matrix τ of whole soft finger can be expressed as

AC

 τ   0 012 

0

 1q4 f j21

01(4 n  2)

 1q3 f

2 j1

0

 3 q4( n1  n2  2)  4 f

2

 2 q4( n1 1)  4

01(4 n  2)

f

1

2 j2

 3 q4( n  n 2) 3 1

2 j3

f

0

2

2 j3

 2 q4( n 1) 3 1

f

2 j2

 014 n  3 

T

(14)

where τ is a 4(N+3)×1 array, and



f j1   q32  q42  , f j 2  q4(2 n1 1)3  q4(2 n1 1) 4 1/2



1/2



, f j 3  q4(2 n1 n2  2)3  q4(2 n1  n2 2)4



1/2

Finally, the dynamic equation of whole soft finger acting by driven moments can be expressed as

10

ACCEPTED MANUSCRIPT

Mq  ΦqT λ  Kq  Q  τ

(15)

2.2. Governing equation of tapping system during contact–impact phase 2.2.1. Governing equation including contact force constraint When the fingertip and rough target surface are in contact–impact phase, the dynamic equation of whole soft finger will change into

Mq  ΦqT λ  Kq  Q  τ  Fc

CR IP T

(16)

where Fc represents contact force vector, which includes the normal and tangential contact forces. Its specific expression is

Fc   0 

01(4 N 8)

Ft

Fn

0

012  

T

(17)

AN US

Here, Fc is a 4(N+3)×1 array. The normal contact force Fn and tangential contact force Ft can be calculated by the LPM (see Fig. 2b) and Coulomb friction law. 2.2.2. Lumped–parameter model

AC

CE

PT

ED

M

For the fingertip contacting with the rough target surface (see Fig. 2a), the contact points C and C' are located at the fingertip and the rough target surface, respectively. The absolute displacement of C and C' are U C  [U CX U CY ]T and UC  [U CX U CY ]T , respectively. The relative displacement between C and C' is uCC  UC  UC  [uX uY ]T . In this paper, the rough target surface is treated as a rigid surface. Hence, it has U CX  U CY  0 , uX  U CX and uY  U CY . Meanwhile, there are equal but opposite contact forces F  [ Ft Fn ]T and F   [ Ft Fn]T at C and C', respectively. In the LPM (see Fig. 2b), a bilinear compliant element in the normal direction and a linear compliant element in the tangential direction are used to describe the deformation of the local contact zone. These elements are connected to a massless particle P which can either stick or slide on the rough target surface. The absolute displacement of the particle P is U P  [U PX U PY ]T , the relative displacement between C' and P is uPC  U P  UC  [uX uY ]T , and the relative displacement between C and P is uCP  UC  U P  [u X uY ]T . Let vt  u X and vn  uY in order to conveniently express the rate of extension of the compliant elements, Vt  uX and Vn  uY in order to conveniently express the rate of relative velocity between C and C', VP  uX in order to conveniently express the tangential velocity of the particle P. The normal compliant element has stiffness k during an initial period of compression and subsequently, a stiffness ke*2 during the restitution period, where e* is energetic coefficient of restitution. k is measured by experiment or calculated by FEM. The formulas to calculate the contact forces are derived by following two cases. 1° For the case of that the normal contact process only experiences single transition between compression and restitution, Shen and Stronge [15] indicate that the normal contact force Fn is related to the small extensions uY of the normal compliant element during the compression period (0 ≤ t ≤ tc1) and the restitution period (tc1 ≤ t ≤ 11

ACCEPTED MANUSCRIPT

tf).

CR IP T

  0  t  tc1  compression  kuY Fn   (18) 2 2  ke u  k e  1 u t  t  t  restitution     * Y * Y c1 c1 f   where tc1 is the transition time from compression to restitution, tf is the separation time. uYc1 is the maximum of the normal compression, which is equal to k 1Fn  tc1  . 2° For the case of that the normal contact processes experience multiple and repeated (i.e. m times) transitions between compression and restitution, a multiple loading–unloading hysteresis loop in the normal direction (see Fig. 4) is developed in this paper. And the corresponding computational formula of Fn is derived as follows, m  2 uYc(i 1)  uYr(i 1)   ku  k e  1 (tr(i -1)  t  tci  i th compression)    Y *  i 2  (19) Fn   m  th ke2u  k  e 2  1 u  u  u  * Yc1  Yci Yr(i 1)   (tci  t  tri  i restitution)  * Y i2  

AN US

where tr(i–1) (i=2, …, m) is the transition time from the (i–1)th restitution period to the ith compression period. tci is the transition time from the ith compression period to the ith restitution period. It should be pointed out that when i is equal to m, trm is the termination time of the frictional impact. uYr(i–1)=–k–1Fn(tr(i–1)) and uYci=–k–1Fn(tci) are the deformation of the normal compliant element at time tr(i–1) and tci, respectively. F t, F n

M

Fn

ED

com-p: compression res-t: restitution

ke*-2

Ft

CE

PT

k

O

kε-2 (1-e 2)(uYc1+uYc2+···+uYcm-uYr1-···-uYr(m-1)) *

u X, u Y

AC

Fig. 4. Force–deformation relations for normal and tangential compliant element during multiple transitions between compression and restitution.

If the particle P slides on the rough target surface, the Coulomb law of dry friction can be used to represent the friction effect. If the particle P is in stick state, Shen and Stronge [15] indicate that Ft is related to the small extensions u X of the tangential compliant element.

 s  Fn Ft   2 k u X

slip stick

(20)

where s  sgn(VP ) represents the motion direction of the particle P. μ is the coefficient of friction between the particle P and the rough target surface. ε–2 is the ratio of tangential to normal stiffness [15]. 12

ACCEPTED MANUSCRIPT

3. Event driven scheme and computational strategy 3.1. Event driven scheme

PT

ED

M

AN US

CR IP T

Because of the frictional constraint, the fingertip will experience multiple stick–slip transitions in the tangential direction. The modelling of stick–slip phenomena belongs to a classical problem of nonlinear dynamics but is still developing and attracting many researchers. To analyze the frictional contact of the mechanical systems, the optional methods [28] can be classified as follows: (a) regularization method; (b) event driven scheme; (c) time stepping method. Awrejcewicz [29–30] apply Hénon’s method to obtain numerical estimations of the stick–slip transitions existing in the Filippov–type discontinuous dynamical systems with dry friction. Kudra and Awrejcewicz [28] simulate a two degree of freedom plane disk, performing one– dimensional translational and rotational motion on the moving belt. The special event–driven scheme together with a numerical simulation algorithm is given. In particular, the transition conditions between the stick and slip are defined. In this paper, the event driven scheme is also used to solve the non–smooth dynamic problem during the tapping event. In the event driven scheme, the governing equation of tapping system will be changed before and after the transition times of stick–slip, compression–restitution and contact–separation. Hence, it is crucial to detect the contact states and the transition times. Here we assume the contact begins at time t=0. The following criteria are given to determine the initial state, the transition time tc from compression to restitution for each compression–restitution loop (see Fig. 4), the transition time tr from restitution to compression, the transition time t01 from stick to slip, and the transition time t10 from slip to stick. (1) Initial tangential state At initial time t=0, the contact is in compression and either initial stick or slip. An upper bound [15] on the angle of incidence for initial stick is

Vt (0) / Vn (0)   2

(21)

AC

CE

(2) Compression or restitution Compression terminates and restitution starts at the transition time tc when the normal relative velocity vanishes, i.e.

vn  t     0 and vn  tc   0

(22)

where δ is an infinite small positive value. Compression terminates and restitution restarts at the transition time tr when the normal relative velocity vanishes.

vn  t     0 and vn  tr   0

(23)

(3) Stick or slip During stick |Ft(t)|≤μFn(t). If the system is in stick state, the transition time from stick to slip t01 can be obtained by finding the time when the ratio of tangential to normal contact force equals the coefficient of friction. In addition, stick–slip transition requires that at time t01, continuation of stick would result in a ratio of forces outside 13

ACCEPTED MANUSCRIPT

the cone of friction.

Ft (t01 )   and Ft (t01 ) Ft (t01 )  0 Fn (t01 )

(24)

Alternatively, if the system is in slip state, the sliding speed of partivle P vanishes at the transition time from slip to stick t10.

VP (t10 )  Vt (t10 )  vt (t10 )  0

(25)

CR IP T

Initial position

Large overall spatial motion

Y

|y|
AN US

|Vt(0)/Vn(0)|<με2

N

Y

Stick compression Eqs (18-20)

Y

vn(t)<0

Stick restitution Eqs (18-20)

M

Y |F (t)/F (t)|<μ t n

Slip compression Eqs (18-20)

N vn(t)<0

N

|Ft(t)/Fn(t)|<μ

Y

Vt(t)-vt(t)=0

AC

CE

PT

N

Y

N

ED

vn(t)<0

N

Y

N

Slip restitution Eqs (18-20)

Vt(t)-vt(t)=0

Y

Y

N Y Recompression

vn(t)<0 N Fn(t)=0

N

Y Separation

Fig. 5. Event driven scheme to determine governing equation.

(4) Separation At time tf–δ, the normal contact force is positive Fn(tf–δ)>0. The contact terminates at time tf when the normal contact force vanishes.

Fn (tf )  0

(26)

14

ACCEPTED MANUSCRIPT

According to the above criteria, the computational flowchart of event driven scheme is shown in Fig. 5. The word ‘slip compression’ in the box means that the motion state of fingertip is slip in the tangential direction as well as compression in the normal direction. Through the event driven scheme, the stick–slip and compression–restitution transition can be determined. Furthermore, the corresponding governing equation can be chosen for the current motion state. 3.2. Computational strategy

CR IP T

After choosing the governing equation, generalize–α method is applied to integrate the governing equation (i.e. Eq. 16). The generalized displacement q and generalized velocity q at time t+Δt are calculated by the following recurrence formulas:

q(t  t )  q(t )  q(t )t   0.5    t 2a(t )  t 2a(t  t )

(27)

q(t  t )  q(t )  t 1    a(t )  ta(t  t )

(28)

AN US

where Δt is the size of time step, β=(1+αf –αm)2/4 and γ=0.5+αf –αm. Vector a should satisfy the following relation:

1   m  a(t  t )   ma(t )  1   f  q(t  t )   f q(t ) ,

a(0)  q(0)

AC

CE

PT

ED

M

where αm=(2ξ–1)/(ξ+1), αf=ξ/(ξ+1). ξ is the spectral radius of algorithm and ξ[0, 1]. ξ=1 will lead to generalized–α method, which is an energy–preserving algorithm, degenerate to be a gradient algorithm. ξ=0, the maximum energy will be dissipated by the algorithm [31]. The generalized–α method will degenerate to be Newmark method when αm=αf=0 or HHT method when αm=0 and αf≠0. In the present study, based on the Arnold’s work [32], the Broyden–Newton [33] approach is introduced to the iteration process of generalized–α method. The corresponding entire computational scheme is illustrated in Fig. 6. In the process, the parameters ˆ  1   m  / 1   f  t 2 and ˆ   / t , which can satisfy q / q  I ˆ and q / q  I ˆ . Here, the variable tol is the integration error tolerance, J is the Jacobian matrix.  Μ ˆ   ΦqT λ    Kq  Q  τ  Fc  ΦqT  G q q  J  T     q λ Φ 0 q   where

 Μq  ΦqT λ  Kq  Q  τ  Fc  G  Φ (q, t )  

15

ACCEPTED MANUSCRIPT

Parameters Ew, Aw, ρw, dw Lw, e* , t ξ, tol, μ, g, B

Initial conditions

q(0)  0 , q(0)  0 , q(0)  0

Start(t=0) M, K(q), Q, τ(q)

q, q, q

q(t  t )  q(t )  q(t )t  (0.5   )t 2a q(t  t )  q(t )  t (1   )a λ(t  t )  0 , an1   f qn   ma  / 1   m 

q(t  t )  q(t  t )  t a q(t  t )  0 dif=1

CR IP T

q(t  t )  q(t  t )  t 2  a

 Mq(t  t )  ΦqT (t  t ) λ  Kq(t  t )  Q  τ (t  t )  Fc ((t  t))  r (t  t )    Φ  q(t  t ), t   

G (t  t )   q(t  t )

λ(t  t ) 

T

AN US

J

r (t )  r (t  t ) , Δ   q λ  J 1r (t  t ) T

q(t  t )  q(t  t )  q(t  t ) , q(t  t )  q(t  t )  ˆq(t  t ) q(t  t )  q(t  t )  ˆq , λ(t  t )  λ(t  t )  λ

M

dif  rn1  tol ?

Y

ED

q(t  t ) , q(t  t ) , q(t  t ) , Fc (t  t )

Y

output

If t
CE

PT

a  a  1   f  / 1   f  q(t  t )

N

Fig .6. Numerical integration procedure based on generalized–α method.

AC

4. Numerical examples In this section, the convergence, accuracy and efficiency of HCM are demonstrated using various tapping systems. The HCM solution and LS–DYNA solution are compared in Section 4.1. Then, based on HCM, the influence of coefficient of friction and structural compliance on the contact forces for the soft finger tapping system without driven moments is analyzed in Section 4.2. In Section 4.3, the movement posture, stick–slip motion state, deflection of phalanx and structural stress of a soft finger tapping system driven by external moments are discussed. In the computation based on HCM, the event driven scheme given in Section 3.1 is adopted to determine the stick–slip and compression–restitution transitions. The contact forces are calculated by LPM. The generalize–α method is used to integrate 16

ACCEPTED MANUSCRIPT

the governing equation of system, and the computational parameters are chosen as ξ=0.8 and tol=1×10–5. Meanwhile, the size of time step is set as t  1ms for the non–impact phase as well as t  0.1μs for the contact–impact phase. All numerical computations are dealt within a Lenovo computer (Ideapad Y700 with Intel(R) core(TM) i7–6700HQ CPU @ 2.60GHz processor, 8GB RAM and Windows 10 operating system). 4.1. Convergence, accuracy and efficiency

CR IP T

The physical parameters of the tapping system are shown in Table 1. Three phalanges are the same. The gravity acceleration is g=9.8m/s2. Unless noted otherwise, H=0.1m as well as μ=0.5 in the following computation. Table 1 Parameters of finger tapping system Ew (GPa)

ρw (kg/m3)

dw (mm)

σw

Lw (mm)

6

k (×10 N/m)

e*

Hard material Soft material

20 1

1048.7 1400

15 15

0.4 0.4

40 40

5.15 0.78

1 1

AN US

Tapping system

Fig. 7 is the contact forces under different N. From Fig. 7, it can be found that the contact forces under N=15 have a good agreement with those under N=60. It shows that the HCM has a good convergence to the density of mesh. Hence, N is chosen as 60 in the following calculation. N=15 N=30 N=60

M

Fn

ED

100

0

-100 0.00

0.03 Time (ms)

Ft

0.06

60

Contact forces (N)

(a)

PT

Contact forces (N)

200

N=15 N=30 N=60

(b)

30

Fn

0

Ft -30 0.0

0.1

0.2

Time (ms)

CE

Fig. 7. Contact forces: (a) hard finger tapping system; (b) soft finger tapping system.

AC

To validate the accuracy and efficiency of the HCM, the HCM solution is compared with LS–DYNA solution. A three–dimensional finite element model (see Fig. 8) is built by LS–DYNA software. The LS–DYNA model’s physical parameters including Young’s modulus, density, Poisson’s ratio, geometric size, gravity acceleration and coefficient of friction are the same as those of HCM. For LS–DYNA model, the solid 164 element is used to mesh it, and the total number of elements is 60261. To simulate the connection constraint of the joint, 4 nodes on the surface of the shaft hole of the phalanges are selected to form a ‘re–joint’ element, where each phalanx provides 2 nodes. The numerical integration method in LS–DYNA is explicit central difference method. The stick–slip transitions and the contact forces are calculated by penalty function method. The factor of time step is set as 0.9. The total integration time is 0.14s, which is the same as those of HCM computation. 17

ACCEPTED MANUSCRIPT

Fixed hinge

Fig. 8. LS–DYNA model.

Fn

100

0

-100 0.00

Ft 0.03

0.06

0.09

70

Contact forces (N)

LS-DYNA solution HCM solution

(a)

(b)

35

0

LS-DYNA solution HCM solution

Fn

AN US

Contact forces (N)

200

CR IP T

Movable hinge

-35 0.0

Ft

0.1

0.2

0.3

Time (ms)

Time (ms)

Fig. 9. Contact forces: (a) hard finger tapping system; (b) soft finger tapping system.

ED

M

Fig. 9 is the comparison of contact forces Fn and Ft between the HCM solution and the LS–DYNA solution. From Fig. 9, it can be found that the HCM solution has a good agreement with the LS–DYNA solution. The normal relative motion of the fingertip for two tapping systems only experience one from compression to restitution transition.

PT

Table 2 Peak value of contact forces

CE

LS–DYNA solution HCM solution Error

Hard finger

Soft finger

Fn (N)

Ft (N)

Fn (N)

Ft (N)

173.104 173.084 0.012%

-86.395 -86.542 0.170%

64.126 64.435 0.48%

-29.683 -32.218 7.80%

AC

Table 2 is the peak value of Fn and Ft. Be in comparison with the LS–DYNA solution, it can be found that the errors of Fn for the hard finger and the soft finger are merely 0.012% and 0.48%, respectively. The errors of Ft for hard finger and soft finger are only 0.17% and 7.80%, respectively. They all show that the accuracy of the HCM solution is high. It also can be found that the soft finger has a good impact resistance. Its peak value of Fn is roughly 30% by those of the hard finger. Table 3 is the comparison of computational efficiency between the HCM and the LS–DYNA model. For the hard finger, the time cost is 21473s for LS–DYNA model while 247s for the HCM. Similarly, for the soft finger, the time cost is 14536s for the LS–DYNA model while 480s for the HCM. All results show that the HCM has a

18

ACCEPTED MANUSCRIPT

much higher efficiency. Table 3 Comparison of computational efficiency Hard finger LS–DYNA solution HCM solution

Soft finger

TC(s)

MO(Mb)

TC(s)

MO(Mb)

21473 247

310.5 584

14536 480

310 582

(TC indicates time costs; MO indicates memory occupancy)

4.2. soft finger tapping system without driven moments

CR IP T

In this section, the influence of μ and structural compliance on the contact forces and structural stress are discussed. The variation of structural compliance is achieved 1/2 through changing Young’s Modulus Ew and slenderness ratio  w  Lw  I w / Aw  . 4.2.1. Effect of μ on contact forces and structural stress

Fn 40

0

Ft

-40

compression

0.1

restitution

40

-40

Ft stick

0.0 0.0 2.5 0.0 0.0

0.1

0.1

0.2

Time (ms) 3

μ=0.3

0.2 μ=0.5

0.2 μ=0.7

Stress (MPa)

PT

CE

0 0.0 2.5

AC

Stress (MPa)

2

compression restitution

slip

0.1

Time (ms) (c)

μ=0.9 μ=1.1 μ=1.3

0

-80 0.0

0.2

ED

0.0

Fn

(b)

Contact forces (N)

μ=0.3 μ=0.5 μ=0.7

(a)

M

Contact forces (N)

80

AN US

The physical parameters of the system shown in Table 1 are adopted in this section. The contact forces and stress are calculated using the HCM.

(d) μ=0.9

0 0.0 3

0.1

0.2 μ=1.1

0 0.0

0.1

3

0.2 μ=1.3

0

0.1

0.2

0.0

0.1

0.2

Time (ms)

Time (ms)

Fig. 10. Contact forces and structural stress at midpoint R of the 3 phalanx under different μ. rd

Figs. 10a and 10b are curves of Fn and Ft, respectively. Fig. 10a shows that the peak value of Fn will decrease as μ increases, on the contrary, the peak value of Ft will increase as μ increases. That’s because there is a large increment of μ, which leads to Ft=μFn still increases even though Fn decreases. Comparing with those of μ=0.3, the peak value of Fn for μ=0.5 and μ=0.7 will decrease 7.2% and 11.6%, respectively. The peak value of Ft for μ=0.5 and μ=0.7 will increase 55.6% and 105.3%, respectively. It 19

ACCEPTED MANUSCRIPT

CR IP T

also shows that it is in gross slip when μ≤0.7. Be different from Fig. 10a, Fig. 10b shows that the peak value of Fn maintain as a constant when μ≥1.1. Meanwhile, there exists a transition from stick to slip. This transition time of stick–slip for μ=0.9 is during the normal compression period. But for μ=1.1 and μ=1.3, the transition time of stick–slip is during the normal restitution period. Figs. 10c and 10d are the normal stress σR of cross section at the midpoint R of the rd 3 phalanx. σR consists of two componients: one is contributed by axial externsion– compression deformation; the other is contributed by benging deformation. The midpoint R is located at the intersection between the cross section and the lower surface of the phalanx. Fig. 10c shows that the peak value of σR will increase as μ increases. Comparing with those of μ=0.3, the peak value of σR for μ=0.5 and μ=0.7 will increase 13.1% and 24.9%, respectively. From Fig. 10d, it can be found that when μ≥0.9, the peak value of σR maintains as a constant even though μ increases. 4.2.2. Effect of structural compliance on contact forces and structural stress

80

Ew=0.1GPa

Fn

Ew=1GPa

M

30

0

-30

Ft

0.0

compression restitution restitution compression

0.2

Contact forces (N)

Ew=0.01GPa

(a)

ED

Contact forces (N)

60

AN US

The structural compliance depends on martial property Ew and geometrical property δw. The variation of δw is achieved through changing dw. The physical parameters ρw, σw, Lw, k and e* illustrated in Table 1 are adopted in the following computation. (b)

Fn

40

0.0 0.0

0.2

0.30 Ew=0.1GPa 0.4 Ew=1GPa

0.1

0.2

(d)

0.2

Stress (MPa)

0 0.0 2.5

0.1

3 Ew=0.01GPa

0.15

δw=10.7

Time (ms)

PT

CE 0.0 0.00 1

AC

Stress (MPa)

(c)

δw=16

Ft

Time (ms)

0.3

δw=32

0

-40 0.0

0.4

compression restitution restitution compression

0 0.00 3.5 0.0 0.0

δw=32 0.05

0.1

2 0 0.0

Time (ms)

0.1

0.10 δw=16 0.2 δw=10.7 0.2

Time (ms)

Fig. 11. Contact forces and structural stress at midpoint R of the 3rd phalanx under different structural compliance: (a) δw=10.7 and μ=0.5; (b) Ew=1GPa and μ=0.5;(c) δw=10.7 and μ=0.5; (d) Ew=1GPa and μ=0.5.

Fig. 11 gives the contact forces and structural stress under different Ew and δw. 20

ACCEPTED MANUSCRIPT

AN US

35

Ew=0.01GPa

0 0.00

0.03

Ew=0.1GPa Ew=1GPa

0.06

0.09

Normal xontact force (N)

M

(a)

ED

70

PT

Normal contact force (N)

CR IP T

From Fig. 11, one can make the following observations:  A large structural compliance will lead to more compression–restitution transitions during the contact–impact period. The fingertip will experience four compression–restitution transitions for Ew=0.01Gpa, two transitions for Ew=0.1GPa as well as one transition for Ew=1GPa. Similarly, the fingertip will experience three compression–restitution transitions for δw=32, two transitions for δw=16 as well as one transition for δw=10.7.  The first peak value of Fn and Ft will decrease as the structual compliance increases. The first peak vaule of Fn for Ew is equal to 0.1Gpa and 1Gpa are roughly 1.5 times and 2.3 times of those for Ew=0.01GPa, respectively. What’s more, the first peak vaule of Fn for δw=16 and 10.7 are roughly 3.2 times and 6.3 times of those for δw=32, respectively.  The maximum of σR will decrease as Ew increases. The maximum of σR is 0.277MPa, 0.983MPa and 2.246MPa for Ew is equal to 0.01GPa, 0.1GPa and 1Gpa, respectively. Similarly, the maximum of σR is 3.170MPa, 3.552MPa and 2.246MPa for δw is equal to 32, 16 and 10.7, respectively. To conveniently plot the loading–unloading hysteresis loop of the normal contact element, e* is chosen as 0.8. Fig. 12 is the force–deformation relation when e*  0.8 . It clearly displays the multiple compression–restitution transitions during the contact– impact period. It also shows that the transitions between compression and restitution increase more as the structural compliance increases. 70

(b)

35 δw=32 δw=16 0 0.00

CE

Deformation (mm)

δw=10.7 0.03

0.06

0.09

Deformation (mm)

Fig. 12. Force–deformation relation of normal compliant element (blue arrow indicates

AC

compression; red arrow indicates restitution): (a) δw=10.7; (b) Ew=1GPa.

4.3. soft finger tapping system with driven moments In this section, the influence of the driven moments on the tapping events is studied. Three cases of driven moments are chosen as: 1) τ1=0.01N·m, τ2=0.008N·m and τ3=30N·m; 2) τ1=1N·m, τ2=0.8N·m and τ3=3000N·m; 3) τ1=1N·m, τ2=0.9N·m and τ3=5500N·m. The driven moments are all in counterclockwise direction. Here, the soft finger has Ew=0.1GPa, ρw=1400kg/m3, dw=15mm, σw=0.4, k=0.78×106N/m, Lw=40mm and e*  1 . 4.3.1 Movement posture during pre–impact, contact–impact and post–impact phases 21

ACCEPTED MANUSCRIPT

Fig. 13 is the moving process of soft finger under different driven moments. In Fig. 13, T indicates the motion time. At the beginning T=0, the fingers are at rest, and all phalanges are located at the horizontal position. From Fig. 13, it can be found that for making the fingertip get to the rough target surface, the system of cases 1, 2 and 3 will cost 0.0912s, 0.01218s and 0.0105s, respectively. For cases 1 and 2, the postures of soft finger are concave when the fingertip is in contact with the rough target surface. Fig 13c shows that posture of the soft finger for case 3 is convex at the impact instance.

T=0.09654s

-0.09

-0.06

con

-0.09

rough target surface

T=0.01308s

-0.03

0.00

-0.12

-0.09

-0.06

-0.03

ct

pa

st-

nta

T=0.0067s

-0.06

im

ct-

pa ct

im

-0.03

-0.09 rough target surface

T=0.01410s

T=0.0105s

T=0.01218s

0.00

-0.12

-0.09

-0.06

-0.03

0.00

X position (m)

X position (m)

X position (m)

initial position t pac -im pre

T=0s

po

ta c

-0.06

act

pac t t -im

e-i

T=0.0091s

(c) case 3

co

t ac mp

-0.03

T=0.0912s

-0.12

0.00

p o st -imp

rough target surface

Y position (m)

-0.09

ct

t

n ta

-0.06

impa

cti

mp ac

-im pre

post-

T=0.0672s

co

Y position (m)

t

pac

-0.03

initial position

CR IP T

case 2 0.00 (b)T=0s

Y position (m)

initial position

T=0s

pr

0.00 (a) case 1

AN US

Fig. 13. The movement posture of the soft finger tapping system in space.

4.3.2. Contact forces and structural stress during contact–impact phase Fig. 14a displays the curve of Fn and Ft under different driven moments. Fig. 14b is the normal stress σR at midpoint R. From Fig. 14, one can make the following observations:

M

1.4

(b)

case 1

ED

200 0 -200

Ft

compression restitution restitution compression

0.3

CE

0.0

0.6

Time (ms)

case 1 case 2 case 3 0.9

Stress (MPa)

Fn

400

PT

Contact forces (N)

(a)

0.0 0.0 12

0.1

0.2

0.3

0 0.00 10

0.05

0.10

0.15

0.4

case 2 0.20 case 3

0 0.0

0.2

0.4

0.6

0.8

Time (ms)

Fig. 14. Contact forces and structural stress under different driven moments.

For case 2, the first peak values of Fn and Ft are roughly 8.5 times of those for case 1. For case 3, the fingertip will experience three compression–restitution transitions. A ‘reverse slip’ phenomenon is found in case 3. Ft experiences multiple transitions from negative to positive during whole contact–impact phase. This phenomenon cannot be found in hard body system but it will occur frequently in soft body system. The fundamental reason is that the soft material has very large tangential contact compliance. This phenomenon should be paid attention on for the contact problem of soft robot. A larger driven moment will lead to a large structural stress. The maximum of σR are 1.296MPa, 10.963MPa and 8.278MPa for cases 1, 2 and 3, respectively.

AC







22

B C D F G H J K L N O P

ACCEPTED MANUSCRIPT

4.3.3. Bending deflection of phalanx Fig. 15 is the bending deflection of the 3rd phalanx during the contact–impact phase. The axial position of fingertip is at –0.04m in Fig.15. Points A and E are the beginning and termination time of impact, respectively. Points B to D illustrated in Figs. 15a and 15c indicate the transition time between the normal compression and restitution.

C 0.2

E 0.4

Time (ms)

C D

D E

A B

E

0.1

0.2

Fn (N)

B C A B neutral layer

-0.04

-0.02

0.00

Axial position (m)

Fn (N)

6 B C

C D

3 A B

0 neutral layer

ED

-0.04

D

B

C

250

A

0 0.0

E

0.3

0.6

0.9

A B C D E

Time (ms)

D E

M

Deflection (mm)

A B C D E

Time (ms)

C D

3

500

-3

D

A

0.00

Axial position (m)

250

0 0.0

neutral layer

-0.02

C B

D E

6

0

-0.04

500

(b) case 2

CR IP T

0 0.0

9

AN US

0.0

25

A

B C

A B C D E

D

Deflection (mm)

0.8

0.4

B

50

Fn (N)

Deflection (mm)

(a) case 1

(c) case 3

-0.02

0.00

Axial position (m)

Fig. 15. Bending deflection of the 3rd phalanx.

PT

AC



In Figs. 15a and 15c, from A to C, the deflection of the 3rd phalanx will continuously increase as time increases. From C to E, the deflection will continuously decrease as time increases. In Fig. 15b, from A to B, the deflection near the fingertip will increase with the increasing of time, but the deflection of the rest will decrease. From B to E, the deflection will continuously increase with the increasing of time too. The occurrence time of the maximum deflection at the fingertip is the termination time of the first restitution period for all cases. And the maximum deflections of cases 2 and 3 are roughly 8.61 times and 7.82 times of those of case 1, respectively.

CE



4.3.4. Stick–slip motion Fig. 16 is the tangential velocity VP of particle P. According to the criteria of motion states, VP=0 shows the fingertip is in stick state, otherwise it is in slip state. When the posture of the soft finger is concave (i.e. cases 1 and 2), the fingertip is in gross slip. However, when the posture is convex (i.e. case 3), the fingertip experiences three 23

ACCEPTED MANUSCRIPT

case 1 0.0 12 (b)

0.1

0.2

0.3

9 0.00

0.4

case 2 0.05

0.10

0.15

0.20

15.0 (c) 7.5 0.0 -7.5 -15.0

case 3

CR IP T

1.0

Tangential velocity (m/s)

1.5 (a)

slip

stick

AN US

Tangential velocity (m/s)

stick–slip transitions. The ‘reverse slip’ phenomenon can be seen much more clearly in Fig. 16c. The basic reason of this phenomenon is that if the slip of particle P is brought to a halt during the contact–impact event, the tangential contact compliance (i.e. the deformation of tangential compliant element) can subsequently reverse the direction of slip[13]. In the first and third slip period, the fingertip slides in the right direction. But it slides in the left direction in the second and fourth slip period. Incorporating with Fig. 14a, it also can be found that the absolute value of the ratio of Ft/Fn is less than μ during the stick period.

0.0

0.3

slip reverse

0.6

0.9

Time (ms)

Time (ms)

Fig. 16. Tangential velocity of particle P.

5. Conclusions

AC

CE

PT

ED

M

In this paper, a HCM is proposed to analyze the frictional impact dynamic responses during the tapping of soft finger. The feasibility of HCM to analyze the tapping event of soft finger has been validated. The applications of HCM and its computational method are demonstrated using various tapping systems. The study found that HCM has higher computational efficiency than three–dimensional LS– DYNA model. In comparison with the LS–DYNA solution, the errors of the HCM solution for the peak value of the contact forces are less than 8.0%. The numerical results show that the soft finger has a good impact resistance due to the large structural compliance. Comparing with the hard finger, the peak value of Fn of soft finger will decrease by 63%. It also shows that the coefficient of friction and structure compliance all have a significant effect on Fn and Ft. When μ≤0.7, the peak value of Fn will decrease as increases μ, and the peak value of Ft will increase as μ increases. The fingertip is in gross slip. However, when μ≥0.9, the peak value of Fn nearly maintains as a constant as μ increases, while the peak value of Ft will still increase. As the structural compliance increases (i.e. Ew decreases or δw increases), Fn will experience multiple compression–restitution transitions. Besides, the peak value of Fn will increase as the driven moments increase. The deflection of the 3rd phalanx gets to the maximum value when the first restitution period terminates. When the posture of the finger is convex, the fingertip experiences multiple repeated stick–slip transitions. Furthermore, the ‘reverse slip’ phenomenon caused by large tangential contact compliance is captured. The advantages of HCM are 1) It has higher efficiency comparing with the 24

ACCEPTED MANUSCRIPT

conventional finite element method; 2) It can describe the large overall motion, large deformation of whole structure and large contact deformation of contact zone at the same time. The proposed method is not hard to be degraded for solving the most developed soft fingers consisting of single phalanx with joint. More importantly, although our investigations only focus on the soft finger, the approach proposed in this paper could also be extended to the walking or climbing of soft robots without difficulty.

Acknowledgements

References

CR IP T

This research was supported by Robotics and Intelligent Machine Lab., NUST and National Natural Science Funds of China (Grant No. 11572157 and 11302107); these supports are gratefully acknowledged.

[1] M.A. Robertson, J. Paik, New soft robots really suck: Vacuum–powered systems empower diverse

AN US

capabilities, Science Robotics, 2 (2017), eaan6357.

[2] S. Terryn, J. Brancart, D. Lefeber, and et al, Self–healing soft pneumatic robots, Science Robotics, 2 (2017), eaan4268.

[3] I. Hussain, G. Salvietti, G. Spagnoletti, et al, A soft supernumerary robotic finger and mobile arm support for grasping compensation and hemiparetic upper limb rehabilitation, Robotics and Autonomous Systems, 93 (2017) 1–12.

[4] S. Helen, Meet the Soft, Cuddly Robots of the Future, Nature, 530 (7588) (2016).

Engineering, 53 (2017) 19–28.

M

[5] J. Zhang, T. Wang, J. Hong, Y. Wang, Review of Soft–bodied Manipulator, Journal of Mechanical

ED

[6] S. Bochereau, B Dzidek, M. Adams, V. Hayward, Characterizing and imaging gross and real finger contacts under dynamic loading, IEEE Transactions on Haptics, 10 (2017) 456–465. [7] A.A. Lat, S. Nefti–Meziani, S. Davis, Design of a variable stiffness soft dexterous gripper, Soft Robotics, 4

PT

(2017) 274–284.

[8] G. G. Rigatos, Model–based and model–free control of flexible–link robots: A comparison between representative methods, Applied Mathematical Modelling, 33 (2009) 3906–3925.

CE

[9] S. J. Cummins, P. W. Cleary, Using distributed contacts in DEM, Applied Mathematical Modelling, 35 (2011) 1904–1914.

[10] Y. Shen, J. Gu, Research on rigid body–spring–particle hybrid model for flexible beam under oblique impact

AC

with friction, Journal of Vibration Engineering, 29 (2016) 1–7. [11] Y. Shen, X. Yin, Analysis of geometric dispersion effect of impact–induced transient waves in composite rod using dynamic substructure method, Applied Mathematical Modelling, 40 (2016) 1972–1988.

[12] K. L. Johnson, The bounce of ‘super ball’, International Journal of Mechanical Engineering, 111 (1983) 57– 63. [13] W. J. Stronge, Impact mechanics, Cambridge University Press, Cambridge, (2000) 112–113. [14] F. Pfeiffer, Energy considerations for frictional impacts, Archive of Applied Mechanics, 80 (2010) 47–56. [15] Y. Shen, W. Stronge, Painlevé paradox during oblique impact with friction, European Journal of Mechanics/A, 30 (2011) 457–467. [16] H. J. Cha, K. C. Koh, B. J. Yi, Stiffness modeling of a soft finger, International Journal of Control Automation and Systems, 12 (2014) 111–117.

25

ACCEPTED MANUSCRIPT

[17] E. Psomopoulou, D. Karashima, Z. Doulgeri, et al, Stable pinching by controlling finger relative orientation of robotic fingers with rolling soft tips, Robotica, 36 (2018) 1–21. [18] Z. Doulgeri, Y. Karayiannidis, Performance analysis of a soft tip robotic finger controlled by a parallel force/position regulator under kinematic uncertainties, IET Control Theory and Applications, 1 (2007) 273– 280. [19] Z. Doulgeri, S. Arimoto, A position/force control for a robot finger with soft tip and uncertain kinematics, Journal of Field Robotics, 19 (2002) 115–131. [20] H. Yuk, S. Lin, C. Ma, et al, Hydraulic hydrogel actuators and robots optically and sonically camouflaged in water, Nature Communications, 8 (2017) 1–12.

CR IP T

[21] F. Chen, W. Xu, H. Zhang and et al, Topology optimized design, fabrication and characterization of a soft cable–driven gripper, IEEE Robotics and Automation Letters, 3 (2018) 2463–2470.

[22] M. A. Devi, G. Udupa, P. Sreedharan, A novel underactuated multi–fingered soft robotic hand for prosthetic application, Robotics and Autonomous Systems, 100 (2018) 267–277.

[23] W. Schiehlen, R. Seifried, P. Eberhard, Elastoplastic phenomena in multibody impact dynamics, Computer Methods in Applied Mechanics and Engineering, 195 (2006) 6874–6890.

Mathematical Modelling, 55 (2018) 616–636.

AN US

[24] L. Zhang, X. Yin, J. Yang, et al, Transient impact response analysis of an elastic–plastic beam, Applied

[25] Y. Shen, Painlevé paradox and dynamic jam of a three–dimensional elastic rod, Archive of Applied Mechanics, 85 (2015) 805–816.

[26] H. Hu, Q Tian, C Liu, Soft machines: challenges to computational dynamics, Procedia IUTAM, 20 (2017) 10– 17.

[27] M. Berzeri, A. A. Shabana, Development of simple models for the elastic forces in the absolute nodal co–

M

ordinate formulation, Journal of Sound and Vibration, 235 (2000) 539–565.

[28] G. Kudra, J. Awrejcewicz, Bifurcational dynamics of a two–dimensional stick–slip system, Differential Equations and Dynamical Systems, 20 (2012) 301–322.

ED

[29] P. Olejnik, J. Awrejcewicz, Application of Hénon method in numerical estimation of the stick–slip transitions existing in Filippov–type discontinuous dynamical systems with dry friction, Nonlinear Dynamics, 73 (2013) 723–736.

PT

[30] P. Olejnik, J. Awrejcewicz, M. Feckan, Modeling, analysis and control of dynamical systems with friction and impacts, World Scientific Publishing, 2018.

CE

[31] J. Chung, G. M. Hulbert, A time integration algorithm for structural dynamics with improved numerical dissipation: The Generalized–α Method, Journal of Applied Mechanics, 60 (1993) 371–375. [32] M. Arnold, O. Brüls, Convergence of the generalized–α scheme for constrained mechanical systems,

AC

Multibody System Dynamics, 18 (2007) 185–202. [33] C. Broyden, A class of methods for solving nonlinear simultaneous equations, Mathematics of Computation, 19(1965) 577-593.

26

ACCEPTED MANUSCRIPT

Appendix Nomenclature nw

Number of discrete elements of the wth phalanx

N

Number of discrete elements of whole soft finger

H

Falling height of soft fingertip Global position vector of an arbitrary point on the sth element

(s)

r

Nodal coordinates vector of the sth element

qe( s )

Generalized acceleration vector of the sth element

S

Shape function of the sth element

S

First derivatives of the shape function S with respect to x

S

Second derivatives of the shape function S with respect to x

T

Kinetic energy of the sth element

M e( s )

Constant mass matrix of the sth element

 (s)

Density of the sth element

U e( s )

Total strain energy of the sth element

U a( s )

Strain energy caused by axial deformation

AN US

(s) e

U

(s) b

Strain energy caused by bending deformation

E

(s)

Young’s Modulus of the sth element

Cross–sectional area of the sth element

I (s)

Moment of inertia of the sth element

 a( s )

Axial strain of the sth element

 (s)

Curvature of the sth element

K a( s )

Axial stiffness matrix of the sth element

K b( s )

Bending stiffness matrix of the sth element

ED

M

A( s )

Total stiffness matrix of the sth element

K e( s )

Generalized external force matrix of the sth element

B Ew

CE

ρw

Acceleration of gravity

PT

Qe( s ) g

CR IP T

qe( s )

Boolean matrix Young’s Modulus of the wth phalanx Density of the wth phalanx

Iw

Moment of inertia of the sth element

Aw

Cross–sectional area of the wth phalanx Length of the wth phalanx

σw

Poisson ratio of the wth phalanx

dw

Diameter of the wth phalanx

δw

Slenderness ratio of the wth phalanx

qp( w)

Generalized coordinates vector of the wth phalanx

qp( w)

Generalized acceleration vector of the wth phalanx

M p( w)

Mass matrix of the wth phalanx

K p( w)

Stiffness matrix of the wth phalanx

Qp( w)

External force vector of the wth phalanx

Φ

Constraint matrix of whole soft finger

AC

Lw

27

ACCEPTED MANUSCRIPT

()q

∂()/∂q

λ

Lagrange multiplier

q

Generalized displacement vector of whole soft finger

q

Generalized velocity vector of whole soft finger

q

Generalized acceleration vector of whole soft finger

M

Mass matrix of whole soft finger

K

Stiffness matrix of whole soft finger

Q

Generalized external force matrix of whole soft finger Generalized moment matrix of the sth element

τ

Generalized moment matrix of whole soft finger

ζj

Rotational angle at node j

Fc

Contact force vector of soft fingertip in LPM

Fn

Normal contact force of soft fingertip in LPM

Ft

Tangential contact force of soft fingertip in LPM

UC

Absolute displacement of point C in HCM

UC'

Absolute displacement of point C' in HCM

UP

Absolute displacement of particle P in HCM

vt

Rate of extension of the tangential compliant element in LPM

vn

Rate of extension of the normal compliant element in LPM

Vt

Tangential rate of relative velocity between C and C'

Vn

Normal rate of relative velocity between C and C'

VP

Tangential velocity of particle P

e*

Energetic coefficient of restitution

k

Stiffness of normal compliant element during compression period

μ

Coefficient of friction

M

ED

J

Jacobian matrix

ξ

spectral radius of generalized–α algorithm

PT

Error tolerance in the generalized–α algorithm

AC

CE

tol

CR IP T

τ

AN US

(s)

28

ACCEPTED MANUSCRIPT

TABLE CAPTIONS Table 1 Parameters of finger tapping system. Table 2 Peak value of contact forces.

FIGURE CAPTIONS Fig. 1. Soft robotic finger taps rough target surface

CR IP T

Table 3 Comparison of computational efficiency

Fig. 2. Hybrid computational model for soft finger tapping system with frictional impact: (a) Flexible

AN US

multi–link; (b) LPM; (c) Discretization of flexible multi–link. Fig. 3. Global positions of nodes on the sth beam element.

Fig. 4. Force–deformation relations for normal and tangential compliant element during multiple transitions between compression and restitution.

M

Fig. 5. Event driven scheme to determine governing equation.

ED

Fig .6. Numerical integration procedure based on generalized–α method.

PT

Fig. 7. Contact forces: (a) hard finger tapping system; (b) soft finger tapping system. Fig. 8. LS–DYNA model.

CE

Fig. 9. Contact forces: (a) hard finger tapping system; (b) soft finger tapping system.

AC

Fig. 10. Contact forces and structural stress at midpoint R of the 3rd phalanx under different μ. Fig. 11. Contact forces and structural stress at midpoint R of the 3rd phalanx under different structural compliance: (a) δw=10.7 and μ=0.5; (b) Ew=1GPa and μ=0.5;(c) δw=10.7 and μ=0.5; (d) Ew=1GPa and μ=0.5. Fig. 12. Force–deformation relation for normal compliant element (blue arrow indicates compression; red arrow indicates restitution): (a) δw=10.7; (b) Ew=1GPa. Fig. 13. The movement posture of the soft finger tapping system in space.

29

ACCEPTED MANUSCRIPT

Fig. 14. Contact forces and structural stress under different driven moments. Fig. 15. Bending deflection of the 3rd phalanx.

AC

CE

PT

ED

M

AN US

CR IP T

Fig. 16. Tangential velocity of particle P.

30