基于SPH/FEM方法的机器人脚在土壤中的沉陷特性研究

夏 婷，何 钢，李晓龙，胡 鹏，陈 新

(河海大学 机电工程学院，江苏 常州 213022)

基于SPH/FEM方法的机器人脚在土壤中的沉陷特性研究

夏 婷，何 钢，李晓龙，胡 鹏，陈 新

(河海大学 机电工程学院，江苏 常州 213022)

基于SPH/FEM方法，将数值模拟方法与贝克经验公式进行了比较；建立了机器人脚与土壤的相互作用模型，进行了机器人脚在土壤中的沉陷研究。结果表明，当机器人脚底部圆弧曲率半径R为16cm时，机器人脚的几何形状更为理想。

SPH/FEM方法；机器人脚；沉陷；数值模拟；软土

随着机器人技术的飞速发展，机器人已经广泛应用于军事、海底和野外空间探测以及事故救灾等人类不能作业的领域[1]，而我国的地理环境复杂，野外工作的机器人将面临着软土等复杂的作业环境，这就很容易发生沉陷量大的问题，不仅影响支撑腿的支撑能力，而且对腿式机器人的平稳性等性能具有重要影响。

目前，不少学者对土壤的沉陷问题进行了研究，方俊等[2]针对用车轮滑转会导致土壤沉陷问题，建立了颗粒-车辆地面离散体动力学模型，对弹性车轮在沙土上因滑转引起的沙土沉陷问题进行了计算。日本学者Shunsuke Komizunai等[3]进行了腿式机器人在松散土壤上的沉陷研究，他们提出了反映脚的沉陷、滑移以及脚与软土间的作用力的接触模型，对双腿机器人动态沉陷和滑移进行了实验模拟和数值验证。Yu等[4]基于土壤的压力-沉陷和剪切模型，建立了轮腿式月球车中腿-月壤相互作用的力学模型，并给出沉陷、脚的结构参数、力学性能和滑移率之间的关系曲线。徐中华等[5]从准静态和动态仿真两方面对利用有限元法分析土壤的大变形问题进行了总结，并指出了解决土壤大变形问题的重要性。

针对机器人脚在软土中易出现沉陷大的问题，本文采用了SPH/FEM耦合的数值模拟方法，建立了底部带有圆弧曲率半径R的机器人脚在土壤上的行走模型，分析了机器人脚在土壤中的沉陷量、土壤的变形情况、应力分布等，分析结果将为复杂环境下机器人脚几何形状设计提供指导。

1 SPH方法简介

SPH(SmoothedParticleHydrodynamics) 方法是近30年来逐步发展起来的一种无网格方法，现在已经被推广应用到土壤切削、流体动力学、空气动力学等力学分析的方面。

SPH方法存在两个核心问题，分别是函数的光滑近似逼近(Smoothedapproximationofthefunctions)和质点的近似逼近(Particleapproximation)。

a.函数的光滑近似逼近。

式中：f(x)为任意空间变量x的函数；D为x的积分区间；x′为x的导数;s为光滑长度，光滑函数W(x-x′,s)又被称作插值核函数(Interpolationkernelfunction)。

b.质点的近似逼近应用在数值计算中，质点近似逼近表达形式为：

式中：ρ为密度；m为质量;N为质点总数。

2 土壤模型

2.1土壤的屈服准则

土壤中最常用的屈服准则为Mohr-Coulomb屈服准则，Mohr-Coulomb屈服准则不仅有不同的S-D效应(StrengthDifferenceEffect)，而且能反映土体的抗压强度。分析时，选取MAT_FHWA_SOIL即MAT147材料模型作为土壤模型。该模型考虑了孔隙水压力效应、塑性硬化和塑性软化等特征，能够使得仿真结果更加准确，适合于土壤沉陷的动态仿真。Abbo等[6]修正了该模型的Mohr-Coulomb面，用来决定剪切强度。

2.2土壤材料参数

土壤材料MAT147的参数有土壤密度、内摩擦角、体积模量和剪切模量等，相关参数的定义见表1。

3 SPH/FEM数值分析方法的验证

3.1贝克原理

美国著名学者贝克(M.G.Bekker)提出了经典的压力-沉陷关系。他认为一个平板在松软土壤上的静沉陷量取决于作用在压板上的载荷大小，由此得出压力P与沉陷量h之间的关系为：

式中：P为接地压力，kPa；kc为土壤的粘聚力变形模量，kN/mn+1；kφ为土壤的摩擦力变形模量，kN/mn+2；n为土壤的变形指数；h为沉陷量，m；b为压板的宽度，也是接触圆盘的半径，m。

压力P与载荷F的关系为：

式中:A为压板与土壤的接触面积。

3.2贝克经验公式与数值模拟的对比

3.2.1模拟分析相关设置

验证贝克经验公式时，土壤为粘性土，即n=0.7，含水量为20%，kc=15.63kN/m0.7+1，kφ=1 415.5kN/m0.7+2，粘聚力为4.2kPa，内摩擦角为20°[7]。

数值模拟时，设置圆盘的半径r=4cm，厚度H=3cm；取土体的大小为20cm×20cm×10cm，土体中间SPH粒子体积为10cm×10cm×5cm。设置圆盘的密度ρ为7.23g/cm3，弹性模量E为1.17MPa，泊松比ν为0.35。设置圆盘采用单元类型3DTet-Solid168，土壤采用单元类型3DSolid164。圆盘采用四面体形式进行自由网格划分，土壤的有限单元部分则采用六面体映射网格划分，土壤中间SPH粒子的个数为60×60×30。圆盘-土壤相互作用模型如图1所示。

模拟时设置计算的时间为26 000μs，给圆盘施加竖直向下的载荷F，大小为 650N。SPH粒子与圆盘之间采用自动点面接触AUTOMATIC_NODES_TO_SURFACE，SPH粒子与周围土壤有限元之间采用固连接触CONTACT_TIED_NODES_TO_SURFACE_OFFSET。

3.2.2仿真结果

根据式(4)，由P得到相应的压力F值，并通过仿真模拟计算得到相应的沉陷量h值。

贝克经验公式计算和SPH数值模拟得到的曲线图如图2所示，二者沉陷量h均随着压力F的增大而增大，即二者的总体变化规律是一致的。在相同力的作用下，二者的沉陷量有一定的误差，数值模拟的沉陷量h值比贝克经验公式计算的h值要小，这是受仿真模拟中的一些参数设置的影响，表明了用SPH/FEM方法来模拟机器人脚在土壤中的沉陷是合理的。

4 机器人脚在土壤中的沉陷特性研究

设置机器人脚的基本模型如图3所示，H=3cm，d=8cm，r=1cm，b=3cm。为了比较机器人脚底部圆弧曲率半径R对土壤应力分布等的影响，取两个相差较大的R值来进行比较分析。经过对比，选取R=8cm和R=16cm。有限元部分土体的大小为20cm×20cm×10cm，土体中间的SPH粒子部分体积大小为10cm×10cm×5cm。模拟分析的其余设置同上。

4.1机器人脚的沉陷量

计算结果显示，曲率半径R=8cm的机器人脚在土壤中的沉陷量为3.308cm，曲率半径R=16cm的机器人脚在土壤中的沉陷量为3.154cm，二者的计算均能达到收敛。从两者沉陷量的差值可以看出当R=16cm时，机器人脚在土壤中的沉陷量更小，即当曲率半径R=16cm时，有利于减少机器人脚在土壤中的沉陷。

4.2土壤的变形分析

机器人脚的刚度较大，而土壤为弹塑性体，所以只有土壤产生形变和应力。当机器人脚底部的曲率半径R取不同值时，土壤的变形是相同的，现对R=16cm时土壤的变形进行分析。此时，土壤的变形可以划分为3个阶段(如图4所示)：

a.直线变形阶段，对应图4曲线上的oa段，接近于线性关系。机器人脚在压力作用下作用于土壤，此时土壤SPH各粒子的剪应力小于粒子的抗剪强度，土壤整体处于弹性状态。当然，土壤中也有少量的压缩变形，这是由于SPH粒子相互挤紧、土体压缩形成的。因此，该阶段也可以称为压密阶段。

b.局部剪切阶段，对应图4曲线上的ab段。在此阶段，机器人脚底部已经开始破坏土体，土壤变形的速率随着时间的增加而逐渐增大。此时土壤中出现了塑性变形区，即土壤SPH粒子中的剪应力达到了其抗剪强度，随着载荷作用时间的变化，塑性变形区从土壤相互接触底面两边缘点开始逐渐扩大，当两个塑性变形区边缘贯通并形成连续滑动面，土壤变形进入下一阶段。

c.稳定阶段，对应图4曲线上的bc段，此时随着时间的推移，土体发生固结，土的抗剪强度逐渐恢复和提高，使机器人脚获得较大的承载力，机器人脚在土壤中的沉陷达到稳定值。

4.3土壤应力分析

机器人脚在垂直力的作用下在土壤中发生沉陷，所以本文只需考虑土壤Z方向的应力分布。机器人脚底部曲率半径R=8cm，当t为5 000μs和26 000μs时,土壤Z方向的应力图分别如图5(a)、6(a)所示。当t=5 000μs时，土壤最先与机器人脚接触的部分受到直接挤压，因此有向下运动的趋势，该部分粒子Z方向应力向下，而此时其余部分的土壤粒子由于受到直接与机器人脚接触部分土壤粒子的挤压，有向上运动的趋势，Z方向的应力向上。当t=26 000μs时，计算达到收敛，机器人脚与土壤的相互作用已经处于稳定状态，此时整体土壤粒子都由于受到挤压而紧密压缩，粒子间的相互作用力较大，均有向上运动的趋势，所以Z方向的应力值都向上。

通过比较图5和图6可知，当机器人脚底部的曲率半径不同时，土壤粒子的Z方向应力分布情况类似，在相同时刻，R=16cm时土壤的等效应力的最大值点数少于R=8cm时土壤等效应力的最大值点数，这说明当机器人脚底部的曲率半径R=16cm时机器人脚的几何形状更理想。

4.4土壤粒子的流动情况分析

当机器人脚底部曲率半径R=8cm，当t为5 000μs和26 000μs时,土壤粒子的流动情况分别如图7(a)、8(a)所示。在初始阶段，土壤在力的作用下向下运动，先与机器人脚接触的土壤粒子受到直接挤压，该部分的土壤粒子向下流动，而周围部分的土壤受到这部分粒子的挤压，有着向左右下方运动的趋势。随着机器人脚进一步的向下运动，更多的土壤粒子受到压缩，中间部分的土壤粒子顺着机器人脚向下运动，而机器人脚周围的土壤粒子则因受到挤压，形成了类似的土拱效应[8]。

通过比较图7和图8可以知道，当机器人脚底部的曲率半径不同时，土壤粒子的速度矢量分布类似。而在相同时刻，R=8cm时土壤粒子速度矢量的最大值略大于R=16cm时土壤粒子速度矢量的最大值。

由以上比较可知，当机器人脚底面曲率半径R=16cm时，相比R=8cm，机器人脚在土壤中的沉陷量和各方向的应力值更小，因此机器人脚底曲率半径R取16cm更理想。

5 结束语

本文基于SPH/FEM方法分析了机器人脚在土壤上行走时的沉陷特性，首先通过比较贝克经验公式和SPH/FEM方法数值模拟时压力与沉陷量间的关系，得出二者曲线是比较吻合的，验证了SPH/FEM方法数值模拟可以用于本课题的研究。然后对底部圆弧曲率为R的机器人脚在土壤中的沉陷特性进行了分析，比较了R=8cm和R=16cm两种情况。结果表明，当R=16cm时，机器人脚在土壤中的沉陷量和Z方向的应力值更小，此时机器人脚的几何形状较为理想。研究成果对设计机器人脚底部的几何形状具有一定的借鉴意义。

[1] 苏红飞.野外机器人导航避障系统设计[D].南京:南京理工大学，2008.

[2] 方俊，闫民.土壤沉陷特性与车轮滑转关系的离散元研究[J]. 科技导报,2010，28(18):23-26.

[3]KomizunaiShunsuke,KonnoAtsushi,AbikoSatoko,etal.Dynamicsimulationofbipedwalkingonloosesoil[J].InternationalJournalofHumanoidRobotics,2012, 9(4):51-57.

[4]YuXiaoliu,FangLei,LiuJinfu.Interactionmechanicalanalysisbetweenthelunarroverwheel-legfootandlunarsoil[C]//ProcediaEngineeringHarbin.Beijing:ElsevierLtd, 2012.

[5] 徐中华，王建华. 有限元分析土壤切削问题的研究进展[J]. 农业机械学报,2005,36(1)：134-137.

[6]LweisBA.ManualforLS-DYNAsoilmaterialmodel147,Mclean,VA[R].Louisville:FederalHighwayAdministrationResearchandDevelopmentTurnaer-FairbankHighwayResearchCenter, 2004.

[7]WongJY.TheoryofGroundVehicles[M].Hoboken:JohnWiley&Sons,Inc, 2001:166-167.

[8] 刘辉，米海珍，文桃,等.土拱效应理论研究现状及其进展[J].低温建筑技术，2011, 11(3):42-45.

Research on the sinkage behavior of the robot's foot walking on soft soil based on SPH/FEM method

XIA Ting, HE Gang, LI Xiaolong, HU Peng, CHEN Xin

(School of Mechanical and Electrical Engineering, Hohai University, Jiangsu Changzhou, 213022, China)

It compares the numerical simulation result with the empirical formula result based on SPH/FEM method, analyzes the sinkage behavior when the robots walk on the soil. The result shows that the structure of the robot's foot is more perfect at curvatureR=16cm.

SPH/FEM method; robot's foot; sinkage; numerical simulation; soft soil

10.3969/j.issn.2095-509X.2015.04.011

2015-03-13

国家自然科学基金资助项目(51375141)

夏婷(1990—)，女，江苏南通人，河海大学机电工程学院硕士研究生，主要研究方向为数字化设计。

TH114

A

2095-509X(2015)04-0044-05