空间机械臂在燃料加注和在轨组装等复杂的空间任务中扮演了至关重要的角色[1-3]．传统多连杆机械臂在非结构化环境下受自由度与几何尺寸的限制，极易与目标发生刚性碰撞[4-6]．目前，具有良好柔顺性和灵活性的柔性机械臂[7-9]受到广泛关注，其能够在非结构化环境中发挥关键作用[10-13]，可以执行柔顺性操作，以缠绕的方式捕获不同形状和尺寸的目标，与航天员之间的交互更加安全[14-15]，在空间任务中展现了广阔的应用前景[16-17]．以象鼻、章鱼臂及哺乳动物的舌头等肌肉性静水骨骼器官为灵感的仿生柔性机械臂，同时具备了操纵灵活性和抓捕稳定性[18-20]，逐渐成为研究热点．柔性机械臂的运动学分析是对其展开设计分析和控制研究的理论基础，然而其材料非线性和运动多自由度加大了动力学分析的难度．传统分段常曲率模型将连续型机械臂离散为一系列相切的圆弧，建模过程得到简化，取得了一定的应用．但模型存在数值奇异性问题，并且没有很好地分析柔性臂的受力情况，无法精确反映柔性机械臂柔软、灵活的特点[21-23]．Cosserat梁理论描述柔性臂中线每个刚性单元的运动状态，其所得到的连续求解结果可以准确地反映柔性臂姿态．研究者对其求解结果进行验证，证实了其模型的准确性[24-26]．然而传统的Cosserat梁理论中运动学方程未知变量过多、求解过程过于复杂，同时材料模量设定方法缺失，降低了模型的精度，阻碍了模型的推广．针对上述问题，本研究提出了一种无差的Cosserat梁运动学方程简化形式和一种基于多项式拟合的模量取值方法．1 构型设计及基础理论1.1　柔性臂构型本研究的对象是由4根柔性线缆驱动的仿生柔性机械臂，其构型如图1所示．柔性臂由硅胶本体和驱动线缆组成，内部没有刚性结构，线缆由舵机驱动，线缆末端与柔性臂末端通过圆盘固定．驱动线缆均匀布置在柔性臂内，其中yci(s)为线缆li到柔性臂中心线的距离，10.13245/j.hust.231256.F001图1柔性机械臂构型yci(s)=a+(b-a)s/L，式中：a为线缆在柔性臂固定端处到中心线的距离；b为线缆末端固定点到末端面中心的距离；L为柔性臂的总长度．因此柔性臂固定端到线缆li某点的距离uci(s)=u(s)+yci(s)，式中u(s)为柔性臂固定端到中心线某点的距离．本研究设计线缆对称分布，关于该柔性臂的设计细节和制作方法见文献[27]．1.2　Cosserat梁建模理论Cosserat梁建模将柔性臂视作一条连续的线，线上每一个点都是刚性单元，变形不仅可以使单元产生平移，而且伴随着旋转运动．为了便于描述，用局部参考坐标系(t(s),n(s),b(s))来描述每个刚性单元，如图2(a)所示．其中t(s)与梁的中心线相切，n(s)和b(s)在横截面内并相互垂直，局部参考坐标系可以反映每一个刚性单元的材料和几何特性．10.13245/j.hust.231256.F002图2Cosserat梁运动描述为了建立刚性单元之间的相互关系，引入刚性单元的自由度描述，如图2(b)所示．忽略横截面内的变形剪切应力，假设k(s)和ξ(s)表示刚性单元相对于b(s)和n(s)轴的曲率，τ(s)为相对于t(s)的挠率，q(s)为沿着柔性臂长度方向的纵向应变，向量ω=(τ,ξ,k)为刚性单元的自由度描述．基于此可得柔性机械臂在三维空间的运动学描述方程[28]：   dt(s)/ds=k(s)[1+q(s)]n(s)-ξ(s)[1+q(s)]b(s);   dn(s)/ds=τ(s)[1+q(s)]b(s)-k(s)[1+q(s)]t(s);    db(s)/ds=ξ(s)[1+q(s)]t(s)-τ(s)[1+q(s)]n(s);du(s)/ds=[1+q(s)]t(s). (1)要想得到柔性臂的形状描述u(s)，须要先求解形变参数q(s)，k(s)，ξ(s)及τ(s)，这须要对柔性臂进行受力分析，如图3所示．10.13245/j.hust.231256.F003图3柔性机械臂受力示意图文献[29-31]对此进行了详细的说明和推导，在此不再赘述．最终得到由线缆驱动的柔性臂的形变微分方程表达式如下    GI+yc2∑i=14Tiτ'=-k ξ yc2(T1-T2+T3-T4)+k yc(T2-T4)+ξ yc(T1-T3);    EJ+yc2∑i=14Tiξ'=-τk yc2(T2+T4)+τyc(T3-T1)+y˙c(T1-T3);    EJ+yc2∑i=14Tik'=τξ yc2(T1+T3)+τyc(T4-T2)+y˙c(T1-T3);EAq=k yc(T1-T3)+ξyc(T4-T2)-∑i=14Ti, (2)式中：G为材料的剪切模量；I为柔性臂相对于坐标轴t(s)的转动惯量；Ti为线缆li施加的拉力；E为材料的弹性模量；J为柔性臂相对于坐标轴b(s)和n(s)的转动惯量；A为材料截面面积．2 运动学方程的简化Cosserat梁模型的求解核心是坐标系的转换，如式(1)所示，求解柔性臂的位置向量u(s)须要先求解关于t(s),n(s),b(s)的9变量微分方程组．由于三个向量是正交关系，因此方程之间并不是相互独立的．本研究利用各个局部参考系与初始局部参考系(t(0),n(0),b(0))之间的几何关系，进行方程组的简化．下面先分析柔性臂的二维平面运动．图4描述的是当柔性臂在XOY平面内运动时，局部参考坐标系间的几何关系，各局部参考坐标系由初始局部参考系绕b(s)逆时针旋转θ(s)角得到．假设初始局部参考系为g(0)=[t(0),n(0),b(0)]T，则有下式恒成立，g(s)=cos θsin θ0-sin θcos θ0001g(0), (3)即t(s)=[cos θ,sin θ,0]T，n(s)=[-sin θ,cos θ，0]T，b(s)=[0,0,1]T．将式(3)代入式(1)可得dθ/ds=k(s)[1+q(s)]. (4)10.13245/j.hust.231256.F004图4二维局部参考坐标系几何关系式(1)所示的微分方程组便可简化成式(4)所示的单变量微分方程．当柔性臂在三维空间运动时，微分方程组的简化原理与上述原理基本一致．由于u(s)仅与t(s)直接相关，因此局部参考坐标系可以看作是由g(0)依次绕b(s),n(s),t(s)逆时针旋转θ(s),φ(s),γ(s)角得到．图5描述了当柔性臂在三维空间内运动时，局部参考坐标系间的几何关系．假设g(0)的旋转顺序为b(s),n(s),t(s)，则有g(s)=Rt⋅Rn⋅Rb⋅g(0)，(5)式中：Rt=1000cos γsin γ0-sin γcos γ;Rn=cos φ0-sin φ010sin φ0cos φ;Rb=cos θsinθ0-sinθcosθ0001.10.13245/j.hust.231256.F005图5三维局部参考坐标系几何关系将式(5)进一步整理可得局部参考坐标系(t(s),n(s),b(s))的θ(s),φ(s),γ(s)表达式如下：t(s)=cos φcos θcos φsin θ-sin φ;n(s)=sin γsin φcos θ-cos γsin θsin γsin φsin θ+cos γcos θsin γcos φ;b(s)=cos γsin φcos θ+sin γsin θcos γsin φsin θ-sin γcos θcos γcos φ. (6)将式(6)代入式(1)可得关于θ(s),φ(s),γ(s)的微分方程组如下：   dφ/ds=ξ(s)[1+q(s)]cos γ-k(s)[1+q(s)]sin γ;   cos φ(dθ/ds)=k(s)[1+q(s)]cos γ+ξ(s)[1+q(s)]sin γ;   cos φ(dγ/ds)-tan γsin φ(dφ/ds)=τ(s)[1+q(s)]cos φ+k(s)[1+q(s)]sin φsec γ. (7)通过微分方程组(7)中第一、三方程先求出φ(s)和γ(s)，然后将结果代入第二方程求解θ(s)．这样的求解顺序及求解思路避免了欧拉角的奇异性．须要注意的是：线缆驱动的柔性臂所产生的挠率τ(s)接近于0，γ(s)接近于0，因而不会导致公式无意义．在求解精度要求不高的情况下，可以认为γ(s)=0［25，27］，式(7)可以进一步化简为二元方程组．3 基于多项式拟合的模量设定方法要想通过式(1)得到软体臂空间形状的整体描述，须要先求解常微分方程组式(2)，得到参数q(s)，k(s)，ξ(s)及τ(s)．材料弹性模量E及剪切模量G会直接影响到式(2)中微分项系数大小，因此其是影响软体臂空间形状求解结果的关键参数，迫切须要建立一套完整、规范的柔性臂模量取值方法．针对上述问题，本研究提出一种基于多项式拟合的模量设定方法．纵向应变q可以反映截面梁单元变形程度，因此弹性模量可以表示为Em(s)=∑i=0mQiqi(s)，(8)式中：m为多项式的最高次数；Qi(Q0，Q1，Q2)为qi对应的多项式系数．利用式(8)将式(2)中关于q的方程转化为Em(s)A(s)q(s)=C11/[C12+Em(s)J(s)]+C12/[C22+Em(s)J(s)]-∑i=14Ti, (9)式中：C11=b2(T1-T3)2；C12=b2(T1+T3)；C21=b2(T4-T2)2；C22=b2(T4+T2)．选取合适的m值并将式(8)代入式(9)，可得到应变q的解析表达式，即可解得Em(s)表达式．对于本研究构型，m=3已经可以满足要求，下面对此情况的弹性模量表达式进行简要推导．将式(9)变形得到如下形式(qm为Em所对应的纵向应变值)，    Z4E33(s)q3(s)+Z3E32(s)q3(s)+Z22E32(s)+Z21E3(s)q3(s)+Z1E3(s)+Z0=0,式中：Z0=-C11C22-C21C12+(T1+T2)C12C22；Z1=(T1+T2)(C12+C22)J(s)-(C11+C12)J(s);Z21=C12C22A(s);   Z22=(T1+T2)J2(s);Z3=(C12+C22)A(s)J(s);   Z4=A(s)J2(s).进一步化简即可得到求解q3的三次方程K0+K1q3(s)+K2q32(s)+K3q33(s)=0，式中：K0=Z0+Z1Q0+Z22Q02;K1=Z1Q1+Z21Q0+2Z22Q0Q1+Z3Q02+Z4Q03;    K2=Z1Q0+Z21Q1+(2Q0Q2+Q12)Z22+2Z3Q0Q1+3Z4Q02Q1;    K3=Z1Q3+Q2Z21+(2Q0Q1+2Q1Q2)Z22+Z3(2Q0Q2+Q12)+           Q3(Q02Q2+2Q0Q12+2Q0Q2+Q12).基于上述计算过程，弹性模量表达式E(s)便可得到．剪切模量G(s)与其相差一个常数系数，其表达式不用另行推导，可以直接通过下式得到，即G(s)=E(s)/[2(1+v)],式中v为材料的泊松比．4 仿真验证为了对本研究提出的模量设定方法的有效性进行验证，开展基于Ansys的有限元分析实验，以得到柔性臂在空间环境中不受重力及其他环境外力下的弯曲形状．4.1　系数拟合实验所采用的设备为MTS Criterion C43.104微机控制电子万能试验机．试样材料为硅橡胶，型号为EcoFlex-0030．对上述实验结果进行系数拟合，拟合结果如下：a. 多项式拟合，当m=0，1，2，3时，Q0=157.710 0，98.970 0，114.220 0，1.113 3 kPa；Q1=-202.550 0，-68.810 0，-1.155 5 kPa；Q2=307.450 0，0.540 0 kPa；Q3= -4.369 9 kPa．b. Mooney-Rivlin本构模型，C10=19.60 kPa，C01=-0.98 kPa．拟合曲线如图6所示．同时，为了进行有限元验证，选取Mooney-Rivlin本构模型材料，将实验数据导入到Origin软件中，利用最小二乘算法对模型参数C10和C01进行拟合，参数测定实验步骤详见文献[31]．10.13245/j.hust.231256.F006图6应力(σ)-应变实验结果4.2　有限元对比分析将柔性臂材料类型设置为超弹性材料，赋予上节中得到的Mooney-Rivlin超弹性本构模型参数．线缆为普通尼龙材料，弹性模量为2.8 GPa，泊松比为0.3．圆盘为铝合金材料，弹性模量为200 GPa，泊松比为0.3．模型具有规则的几何形状，划分为六面体单元．柔性臂主体网格尺寸为5 mm，共生成2.614 3×104个节点，5 490个单元．模型网格划分情况见图7．10.13245/j.hust.231256.F007图7模型网格划分情况基于上述参数开展柔性臂变形实验分析，并将有限元仿真结果与模型计算结果进行对比．为了更好地验证模型的准确性，取Em(m=0,1,2,3)的仿真结果进行对比分析．以下假设ui(i=0,1,2,3)为当模量取值为Ei(i=0,1,2,3)和Gi时柔性臂的姿态曲线，un为有限元仿真结果．以下假设工况1为T2=3 N，工况2为T2=8 N．图8(a)给出了当线缆拉力取值较小时柔性臂的弯曲状态，曲线几乎重合．图8(b)给出了当拉力取值较大时柔性臂的弯曲状态，此时材料的超弹特性显现，随着应变q的最高保留次数的增加，模型误差减小，进一步证明了方法的有效性和准确性．10.13245/j.hust.231256.F008图8不同工况下柔性臂弯曲形状4.3　模型求解速度对比分析由于本研究提出的运动学方程简化方程的求解结果与原式相同，因此本研究不再给出具体仿真图形．下面通过对比算法简化前后的运行速度，评价简化算法在加快求解速度方面所起作用．本研究使用的算法运行平台为Lenovo Legion Y7000 2020笔记本电脑(Intel i7-10750HQ，2.6 GHz CPU，16 GiB RAM)，算法的运行环境为Matlab R2019a．两种运动学方程在算法编写时保持底层逻辑相同，为了保证结果的可靠性，每个实验重复20次，最后结果取平均值．表1给出了不同工况下柔性臂形状原方程求解时间(t0)和简化方程求解时间(t)的对比结果，提升百分比p=[(t-t0)/t0]×100%．其中，工况1和工况2与上节中设定的工况相同，是对柔性臂在二维平面内的变形进行分析．工况3为T1=3 N，T2=3 N，工况4为T1=5 N，T2=8 N，是对柔性臂在三维空面内的变形进行分析．表1给出的实验结果表明：二维和三维的运动学方程简化形式将模型的求解时间分别减少了约26%和18%，验证了简化方程在提升模型求解效率上的有效性．10.13245/j.hust.231256.T001表1速度实验结果对比工况t0/st/sp/%10.562 20.411 126.8820.553 80.411 825.6430.870 20.709 218.5040.872 00.720 817.345 结语Cosserat梁理论为柔性臂静力学分析提供了连续精确解，针对传统Cosserat梁模型运动学方程未知变量冗余、求解过程过于复杂的问题，本研究基于方向余弦阵，提出了一种无差的运动学方程简化形式，提高Cosserat梁模型的求解速度．同时针对材料模量对于模型的最终求解结果影响较大，而目前没有一个准确且通用的模量设定方法这一问题，提出基于多项式拟合的弹性模量设定方法，可以根据精度要求，调整应变保留项的最高次数；也可以根据不同超弹性材料和变形情况，选取不同的拟合函数，有限元仿真对比实验证明了模量设定方法的有效性．