面向软体多指手的指尖接触力学建模①

陈凌峰 张亚琪 马小龙 蔡世波 鲍官军

(浙江工业大学机械工程学院 杭州 310023)

0 引 言

近年来，软体机器人成为热门的研究方向[1]，使用软材料制成的软体机械手也引起了国内外研究学者的广泛关注和研究[2-6]。以软材料代替刚性材料为解决刚性机械手环境适应性差、交互安全性差等问题提供了新的思路。从生物学角度来看，各种软组织是生物体中不可或缺的组成部分。柔软的组织使生物具有极强的环境适应能力，而且在运动时具有极好的动态特性[7,8]。于是研究学者从仿生学得到灵感，使用软材料来设计和制造软体手，充分发挥软材料固有的柔顺性以及粘弹性和迟滞性等特性，使之具有与生物体相媲美的适应性和灵活性以及交互安全性。因此，软体机械手所具有的复杂环境适应性、交互安全性等性能优势将使其在社会生产和生活的诸多方面发挥重要作用。

软体多指手在本体结构上，尤其是材料的柔性上和人手类似，但是要使软体多指手像人手一样灵活，还需要相应的控制理论和方法。目前的软体多指手研究多数停留在开环抓取研究，仅文献[9]通过位姿解算实现了对目标物体表面形貌的间接感知。传统的机器人多指灵巧手以刚性的关节型机械结构为主，已经形成了较为成熟的理论和方法，例如目标物体聚类分析、抓持位形规划、抓持触点规划、抓持力分配、手指关节力分配、物体操控等。这些理论方法中，目标物体聚类分析、抓持位形规划、抓持触点规划不牵涉到具体的多指手结构和参数，同样适用于软体机械手。然而，由于结构原理的差异，基于传统刚性关节型机器人理论的刚性机械手的抓持力控制和目标物体操控方法并不适用于软体结构的多指灵巧手。文献分析表明，针对软体多指手的目标物体抓持过程交互作用机理、抓持规划理论等方面的研究尚未开展。

手指与目标物体的作用机理和数学模型是进行多指灵巧手抓持规划的力学基础，软体多指手由于其特殊的手指材料和结构，不能直接借用刚性手指的无摩擦点接触或有摩擦点接触模型[10]，为数不多的软指接触模型[11,12]也只是在有摩擦点接触模型的基础上增加了摩擦力矩项，依然是针对刚性结构的多指灵巧手。所以，研究软指尖接触力模型等针对软体多指手与目标物体的作用机理，并建立数学模型对软体手的控制和抓持规划具有重要意义。

软指接触力在抓持稳定性以及安全性和操纵过程的处理中起着重要作用。抓取力的估计需要知道接触特性，包括法向接触力和接触面积之间的关系，以及手指-物体界面处的压力分布情况。Xydas等人[13,14]开发了一个有限元模型，对半球形软指尖进行了有限元分析，并与实验结果进行了比较，以进一步了解软指接触行为。基于非线性弹性有限元建模，他们研究了具有非线性大变形的软指接触。有限元方法可以观察诸如手指接触区域内的接触压力和非线性的力-位移关系等的参数。Kao等人[15]实验证明由变形引起的弹性力与指尖的位移呈幂律关系，并且认为他们的理论包含了赫兹接触。然而，这些研究没有区分软指尖的材料非线性和由指尖的半球形引起的几何非线性，并且只定义了一个同时包含2个非线性效果的参数。因此，难以找出基于模型的模拟结果与实际实验结果之间的差异原因。由于他们提出的模型的复杂性，这些研究不适用于软指操纵系统的运动方程分析。虽然有限元分析可以在软指尖上得到应力分布和弹力，但这些模拟结果取决于所选的单元和网格划分。

本文对软指尖和目标物体的接触作了简化和抽象，建立软指尖和目标物体正向接触的模型。提出了软指尖的微分表示方法，建立了软指尖接触的力学模型，揭示出软指尖与目标物体的接触机理。基于该模型，研究了软指接触力以及接触摩擦力矩的分布情况。在此基础上更进一步研究了任意方向接触的情况，提出了软指尖静态弹性模型。最后对软指尖接触理论进行了有限元分析和实验验证，通过对比分析，表明理论模型的有效性。

1 软指尖正向接触建模

1.1 软指尖正向接触几何模型

本文以半球形来表示软体手指的软指尖，当软指尖和刚性物体接触时，会存在一定的变形。多指手灵巧手的软指尖在接触目标物体时，典型的正向接触如图1所示。

图1 指尖-物体正向接触

图1展示了软指尖沿法向接触刚性目标物体的接触状态。假设软指尖的半径为a，则圆表面上的点满足如下关系：

x2+z2=a2

(1)

其中a为半球半径。指尖表面C点的x坐标在平面中表示为

Cx=(a-d0)tanθ

(2)

其中，d0表示软指尖的变形量，θ表示接触角，如图1所示。联立式(1)和式(2)，C点的z坐标可以表示为

(3)

则CB在z方向的长度表示为

CB=Cz-(a-d0)

(4)

CB的长度在物理上表示软指尖表面上的一点在接触变形时的变形量。最大接触角θmax可以表示为

(5)

通常情况下，软指尖接触面的接触面积随指尖变形量的增大而增大。当变形量为d0的时候，如图2所示，在接触面B点处的半径为

r=(a-d0)tanθ

(6)

对上式两边求微分可以得到半径的微分dr如式(7)所示。

(7)

于是得到因变形而产生的接触面积的微分可以表示如下：

(8)

图2 接触面的微分

1.2 软指尖的微分表示及接触力分析

对于如图1所示的半球形指尖和目标物体的接触模型中，指尖和目标物体形成的接触面形状沿着软指尖的法线对称。为了得到合适的半球形指尖接触力学模型，提出2个关于材料特性的假设。

(1)忽略材料的不可压缩性。

(2)构成软指尖的材料杨氏模量是常数。

使用微积分的思想将软指尖视为无限个柱状微元的组合，每一个柱状微元的轴向沿着z轴方向。图3对于每一个这样的微元进行了力学分析，图4画出了3个这样的柱状微元，微元受到的应力表示为

(9)

其中，F是面积为S的微元受到的接触力的大小。

图3 指尖微元

图4 指尖的微分示意图

当一个小的压缩接触力施加到指尖微元，微元的长度会被压缩，压缩应变ε定义为软指微元长度的微小变化，如式(10)所示。

(10)

其中ε是无量纲的，在受到压缩力的情况下，ε是负值。在机器人抓持任务中，多是属于这种情况。根据软指尖的材料特性和微元的受力分析，使用杨氏模量，微元的受力关系如式(11)所示。

(11)

另一方面，微元受到的压力表示为

(12)

联立式(11)和式(12)，微元上的受力情况最终可以表示为

(13)

从式(13)中可以发现软指微元表面的接触力取决于压力分布、软指微元的长度及其变化量、接触面积和软指尖的杨氏模量。

将软指尖视为一组有无限个软指微元组成的模型，如图4所示，软指尖的变形可以看成是由各个微元变形的组合。由于指尖的力由压力产生，所以所有的微元的压力分布就决定了指尖总的接触力。接触力分布的表达式为

(14)

根据式(14)可以绘制出软指尖的压力分布图，图5(a)、(b)和(c)展示了根据变形量d0分别为1 mm、2 mm和3 mm时的理论接触压力分布情况，从接触力分布图可以看出接触中央的接触力最大，越往接触的边缘接触力越小，直至接触边缘的接触力为零。这种现象很容易从软指尖的微分模型中得到解释，由于接触中心微元的应变最大，相应产生的接触压力也最大，而接触边缘的微元没有压缩变形，故产生的接触压力为零。

(a) d0=1 mm

(b) d0=2 mm

(c) d0=3 mm

对软指尖接触力分布进行积分，可以求得接触面总的接触力Fc。

(15)

其中，E表示指尖的杨氏模量，可以通过实验的方法得到。根据软指尖的变形情况，长度参数L的表达式为

(16)

ΔL的表达式为

(17)

从式(15)可以看出，接触力Fc只是软指尖的变形量d0的函数，并且指尖接触力和变形量d0的平方成正比。

根据软指尖接触力，容易求得软指尖接触所能提供的最大摩擦力Ff为

Ff=Fc·μ

(18)

接触面上各点最大摩擦力对接触面的中心点取矩得到接触摩擦力矩的分布M。

(19)

根据式(19)可以画出如图6所示的接触摩擦力矩分布图。由于是对接触面的中心点取矩，故接触面中心的摩擦力矩为0。在接触面的边缘，由于接触力趋于0，故所能提供的摩擦力为0，因此摩擦力矩也为0。因为软指尖接触可以提供接触摩擦力

图6 接触摩擦力矩分布图

矩，所以当物体受到一个绕轴线的外力矩扰动时，软指尖会产生一个反方向的摩擦力矩与之平衡，从而使软指尖对目标物体的抓持更加稳定。这也是软指尖接触与点接触或摩擦点接触的区别之一。

2 软指尖静态接触弹性模型

对软指尖正向接触模型作推广，更进一步研究软指尖和目标任意方向接触的情况。一般化的指尖-物体接触模型如图7所示，和正向接触模型相比，该模型引入了接触的方向角θp和方位角Φ。为了建立软指尖的弹性力方程，将半球形软指尖视为一组无限数量的虚拟线弹簧组合而成，图7展示了一个这样的线弹簧PR。为了简化方程的推导过程，同样以材料的不可压缩性假设和杨氏模量为常数的假设为前提进行分析计算。

图7 一般化的指尖-物体接触模型

2.1 软指尖刚度

根据接触面的方程xsinθp+zcosθp=a-d，弹性力dF可以表示为

dF=k·PQ

(20)

如图8所示，设k0为弹簧的刚度系数，dS0为截面积，L0是用于计算杨氏模量的柱形样本的原始长度，E是软指尖的杨氏模量。根据线性材料的特性，以下式子成立：

σ=Eε

(21)

从而得到：

(22)

图8 软指尖内的弹簧微元

进一步得到：

(23)

其中，F表示柱形样本受到的压力，δx是压缩的位移。由于假设了杨氏模量是一个不变的物理常数，所以以下的方程是成立的。

(24)

从而k可以表示为

(25)

设K是如图7所示软指尖整个变形部分的刚度，K可以表示为

(26)

其中，

(27)

(28)

其中，ell表示如图7所示的椭圆区域。借助于Matlab对式(26)做一个数值积分，发现不同方向角下的指尖刚度是一个常数，这表明指尖刚度K与目标物体的方向角无关。

如图9所示，设Σ′是从Σ坐标系转换到O′的坐标系，Σ″是z′轴转过θp角度的圆柱坐标系。设r是以点C为圆心的接触圆的半径，φ是绕着z， z′和z″轴的公共旋转角。Σ′坐标系中的(x′， y′)和Σ″坐标系中的(r， φ)的关系可以表示为

x′=rcosφcosθp

(29)

y′=rsinφ

(30)

由于(x， y)和(x′， y′)的关系可以表示为

x=x′+(a-d)sinθp

(31)

y=y′

(32)

所以，Σ坐标系和Σ″坐标系的变量转换可以表示为

x=rcosφcosθp+(a-d)sinθp

(33)

y=rsinφ

(34)

同时，如图9所示的软指尖底面的椭圆区域根据以上的转换规则可以转换为圆形区域，即(r, φ)的积分区域范围分别是[0,ac]和[0, 2π]。

图9 各坐标系的变换关系

于是，令x=rcosφcosθp+(a-d)sinθp，y=rsinφ，则式(26)可以变换为

(35)

由于K是一个和θp无关的常数，因此令式(35)中的θp=0得到：

(36)

2.2 软指尖弹性力

根据式(20)、式(25)和几何关系QT=PQcosθp(如图9所示)，弹性力F可以写成：

(37)

与K的推导过程相类似，对式(37)作变量转换，式(37)变成：

(38)

其中，

(39)

在式(38)中，B(r, φ)对应于式(35)括号中的被积函数，即：

(40)

与式(36)的推导方法类似，最终计算得到弹性力F的表达式为

(41)

3 软指尖接触有限元仿真

有限元法求解的原理是把连续的求解区域分割成有限个单元，并在每个单元上指定有限个结点，假设一个简单的函数近似地表示其位移分布规律，再利用弹塑性理论中的变分原理或其他方法，建立单元结点的力和位移之间的力学特性关系，得到一组以结点位移为未知量的代数方程组，从而求解结点的位移分量。进而利用插值函数确定单元集合体上的场函数。由位移求出应变，由应变求出应力。本文使用Abaqus2017来求解软指尖与目标物体的接触问题，通过有限元仿真结果和上述理论分析结果的对比来验证理论结果的正确性。

Abaqus有限元分析过程可分为建模阶段、计算阶段和后处理阶段。建模阶段是根据结构实际形状和实际工况条件建立有限元分析的计算模型，从而为有限元数值计算提供必要的输入数据。有限元建模的中心任务是结构离散，即划分网格。计算阶段的任务是完成有限元方法有关的数值计算，这部分工作由有限元分析软件控制并在计算机上自动完成。后处理是对计算输出的结果进行必要的处理，并按一定方式显示出来，以便对问题进行分析。

根据图7的理论分析模型建立的有限元分析几何模型如图10所示，对于软指尖与目标物体正向接触的情况，采用1/4球体进行建模。1/4球体代表软指尖，目标物体通过一条直线来模拟。对于非正向接触的情况采用半球模型进行建模，改变直线的倾角以实现不同接触角的加载仿真。在仿真中设置软指尖的半径R为20，分别对接触角θ等于0 °、10 °、20 °、30 °的情况进行仿真。

图10 仿真几何模型示意图

划分网格是建立有限元模型的一个重要环节，有限元网格数量的多少和质量的好坏直接影响到计算结果的精度和计算规模的大小。一般来讲，网格数量增加，计算精度会有所提高，但同时计算规模也会增加，所以在确定网格数量时应该权衡这2个参数。网格较少时，增加网格数量可以显著提高计算精度，而计算时间不会有很大的增加。所以应注意增加网格数量后的经济性[16]。为了更加精确地捕捉模型在变形中微小的变化，对模型中接触区域的网格进行细化。如图11所示，将指尖的接触模型分为3个具有不同网格密度的区域，兼顾了计算的精度和计算的时间。由于接触发生在软指尖的表层，为减少计算时间，同时减小计算误差，软指尖外围的区域具有最大的网格密度。其中，区域I为半球的外围区域，径向范围约为半径的1/4，区域II为半球的中层区域，其径向范围约为半径的1/4，区域III为半球的中心区域，半径为半球半径的1/2。半球体外围具有最大的网格密度，并通过反复地调整网格密度和分析步长使得模型收敛，最终模型由CAX4R四边形单元构成，接触区域最小网格大小为0.2。

图11 软指尖的有限元接触模型网格

软指尖的材料为硅胶，属于超弹性材料，对于这种类型的材料，Abaqus提供了多种本构模型供选择，包括Mooney-Rivlin、Neo Hooke、Yeoh、OGden等。对于较大变形的情况，使用OGden本构模型可以很好地模拟，故在仿真中选用OGden作为本构模型，并输入相应的实验数据定义材料的属性。

在有限元仿真中，边界条件的确定至关重要，边界条件直接影响仿真的精度，甚至会影响收敛性。该模型的边界条件如图12所示。半球底边上的节点在各个方向都是固定的，该边界条件是有效的，原因是软指尖一般固定在手指上不会发生移动。在对称面上的节点进行x方向的约束，由于对称性约束，

图12 模型的边界条件和载荷

对称面上的节点只能在y方向上移动。对于以直线代替的刚性目标物体，要求它只能沿y轴方向平移，所以除了对称性约束，即x方向上的约束之外，还要约束其面内的旋转。

在接触问题的模拟中一般有2种施加载荷的方法，给刚性平板施加位移，然后计算产生的接触载荷；给刚性平板施加压力，然后计算得到位移。在本文的接触模型中，采取的是第2种方法，通过施加不同的压力载荷，来求解相应的接触参数，因为施加压力可以极大地降低所建立模型的收敛难度。并且在对模型进行分析时，如果直接加载相应的压力载荷，整个接触分析可能不会收敛，所以可以通过先施加一个较小的压力载荷，使接触平稳，随后再施加实际的压力载荷。这样虽然会增加分析步的数目，但可以使模型有效收敛，从而减小了计算时间。

建立完指尖接触的有限元模型之后，创建作业并提交分析运行，等待计算完成即可得到有限元仿真的结果。软指尖在一定压力载荷下的变形图如图13所示，从软指尖的接触变形图可以看出，仿真结果的变形情况和实际情况非常相似。

图13 软指尖变形的位移云图

在得到有限元仿真的结果后，需要从结果数据文件中提取需要的数据以便于和理论的结论进行对比。在Abaqus中不能直接获取压力载荷关于软指尖变形量之间的关系数据和曲线，因此，需要对仿真的结果做进一步的后处理。在Abaqus中可以通过Visualization模块中的工具获取节点的位移和时间的关系数据曲线以及压力载荷和时间的关系数据。而借助Visualization模块自带的Combine(X，X)可以将以上两者综合起来得到压力载荷和节点位移之间的关系数据和曲线。软指尖的变形量是通过刚性平板的y方向的位移计算得到的，平板y方向的位移乘上倾角的余弦值即得到指尖的变形量。对于不同的压力载荷，软指尖的变形量不同，通过提取不同压力载荷下的刚性平板位移来计算得到软指尖的变形量。

通过以上的处理得到了在θ=0 °时，软指尖的变形量和压力载荷的关系曲线和数据。通过Matlab将仿真获得的变形量和压力之间的关系曲线和理论分析得到的曲线画在同一幅图上，以便于对比两者之间的关系，如图14(a)所示，图中的虚线表示理论曲线，实线表示有限元仿真曲线。

在本次仿真中还分析了不同接触角情况下的指尖变形量和压力之间的关系，根据图10所示，分别建立接触角为10 °、20 °、30 °的有限元仿真几何模型。和以上平压接触的处理方式类似，通过建模、计算、后处理得到了接触角为10 °、20 °、30 °情况下指尖变形量和接触压力的关系曲线和理论曲线的对比图，分别如图14(b)、(c)、(d)所示。

从4幅理论和仿真的曲线对比图来看，每幅图理论曲线和仿真曲线的关系基本类似。在软指尖的变形量小于5 mm时，对于接触角等于0 °、10 °、20 °、30 °，有限元曲线和理论曲线基本吻合。但是当指尖的变形量超过5 mm时，有限元曲线和理论产生了偏离，在相同的变形量时，仿真的接触压力要大于理论计算得到的接触压力，而且随着变形量的增大，两者的差距越来越大。这是因为在理论分析的时候忽略了材料的不可压缩性，当软指尖的变形量不大时，这种不可压缩性对指尖的变形几乎无影响，但是一旦软指的变形量超过一定的范围，这种材料不可压缩性产生的效应就变得十分明显。这种效应造成的结果是需要使用更大的力才能使软指尖产生少量的压缩变形，这也是当指尖变形量大于5 mm时接触压力快速增加的原因。所以从仿真结果和理论结果的对比分析发现，理论有一定的适用范围，在变形量不是较大以至于不可忽略材料的不可压缩性时，理论都可以适用，而且理论可以比较准确地描述指尖接触力和变形量之间的关系。

(a) 接触角θ=0 °

(b) 接触角θ=10 °

(c) 接触角θ=20 °

(d) 接触角θ=30 °

图14 不同接触角下压力于变形量的关系对比图

4 实验验证

压缩样品采用3D打印的模具制作，如图15所示。软指尖样品的直径分别为20 mm、30 mm、40 mm，弹性模量分别为1 MPa和0.3 MPa。

软指尖的压缩实验在如图16(a)所示的REGER 3010万能试验机上完成，该试验机适用于对橡胶、塑料、金属等材料进行拉伸、压缩等力学性能测试。实验系统具有开放的数据结构，允许用户随机调用结果参数和过程数据。实验数据可方便导入Excel表格，方便用户后期处理。

在进行完标定等准备工作之后，对样品进行压缩实验，如图16(b)所示。在实验中，按照位移的方式加载，根据硅胶材料的性质设置了压缩速率为2 mm/min，根据不同指尖的大小、弹性模量设置合适的加载量以免超过力传感器的量程。对于每一种软指尖进行了3次重复实验。测试完成后，将每组实验的数据导出到Excel表格。通过Matlab画出压力和变形量的关系曲线，如图17和图18所示。

图15 软指尖样品

(a) REGER 3010万能试验机

(b) 样品压缩实验

图17画出了弹性模量为1 MPa，直径分别为20 mm、30 mm、40 mm的软指尖压力和变形量之间的关系，其中的虚线表示理论仿真曲线，从理论仿真结果来看，指尖的接触压力和变形量有关而和指尖的直径无关。根据实验曲线和理论曲线的对比可以发现，直径为40 mm的软指尖实验曲线和理论曲线最为吻合，直径为30 mm的软指尖次之，直径为30 mm的软指尖只有在变形量小于3 mm时与理论曲线重合，当变形量超过3 mm时，实验曲线和理论产生了偏离，在相同的变形量时，实验的接触压力要大于理论计算得到的接触压力，而且随着变形量的增大，两者的差距越来越大。这是因为在理论分析的时候忽略了材料的不可压缩性，当软指尖的变形量不大时，这种不可压缩性对指尖的变形几乎无影响，但是一旦软指的变形量超过一定的百分比，这种材料不可压缩性产生的效应就变得较为明显。这种效应造成的结果是需要使用更大的力才能使软指尖产生少量的压缩变形，这也是当指尖变形量超过一定的百分比时接触压力快速增加的原因。所以从实验结果和理论结果的对比分析发现，理论有一定的适用范围，在变形量不是较大以至于不可忽略材料的不可压缩性时，理论都可以适用，而且理论可以比较准确地描述指尖接触力和变形量之间的关系。

图17 同一弹性模量不同直径的指尖变形量与压力的关系曲线

图18画出了直径为40 mm，弹性模量分别为1 MPa和0.3 MPa的软指尖压力和变形量之间的关系，其中的虚线表示理论仿真曲线。对比不同弹性模量的实验曲线，在相同的变形量时，弹性模量大的指尖所产生的接触压力大于弹性模量小的指尖，这一点符合理论的结果。通过对比实验和理论曲线，无论是弹性模量为1 MPa还是0.3 MPa，在变形量较小(小于30%)时都和理论曲线吻合，在变形量较大时发生偏离，原因与上述的分析相同。此外在软指尖弹性模量较小时，实验曲线和理论曲线吻合的范围更大。

图18 不同弹性模量指尖变形量与压力的关系曲线

5 结 论

本文以软指尖为研究对象，建立了软指尖和目标物体正向接触的模型。提出了软指尖的微分表示方法，有效地解决了非线性接触建模的难点，建立了软指接触的力学模型。探究了软指尖接触机理，研究了软指接触力及接触摩擦力矩的分布情况。并且提出了一般化的软指尖静态接触弹性模型，推导出了软指尖刚度和弹性力的表达式。此外，进行了软指尖接触的有限元仿真及验证实验，结果表明，本文提出的理论模型在软指尖变形量不超过30%时与实际情况吻合良好，但是变形量过大时会产生偏差。这是由于理论分析忽略了材料的不可压缩性，在实际情况下，当软材料的变形量过大时，材料的不可压缩性会变得明显。