探月工程是人类探索宇宙的首要目标，也是迈向深空的第一步．为了深入开展月球科学研究、勘查和利用月球资源，中国已经启动无人月球科研站的立项工作，计划在世界上首次实施无人月球科研站建设的重大工程任务[1-3]．无人月球科研站是由多个协同工作的着陆器、巡视器、功能舱、特种机器人、能源站等月基装备组成的智能集群系统，其建设涉及月基装备的转运、对接、作业等一系列重要活动[4-7]．因此，月基装备的安全转运是无人月球科研站建设过程中的关键科学问题．然而，高内摩擦角、低黏聚力月壤与月基装备足垫的耦合作用机理尚不明确，导致月基装备的转运稳定性难以保证[8-9]．国内外许多学者从理论建模、仿真分析和试验验证等方面研究了足垫与月壤的耦合作用[10-11]．Yin等[12]对足垫与星壤风化层的法向接触力进行了理论和离散元仿真研究，并分析了星壤参数对法向接触力的影响．Witte等[13]进行了足垫-月壤接触过程的多体动力学数值仿真，提出了月基装备着陆点选择方法．钟世英等[14]通过界面直剪试验研究了模拟月壤与足垫材料的接触界面特性，分析了接触面的变形和强度特性．这些研究多聚焦于足垫-月壤的竖向刺入力学，对月基装备行走过程中足垫的水平滑移研究还较少[15-17]；同时，现有的研究多依赖于经验公式或试验基础，考虑月壤力学性质的足垫-月壤接触力学理论分析也较少．本研究提出了一种基于月壤极限承载力理论的月基装备足垫-月壤水平拖曳理论模型．考虑月壤的力学特性，采用对数螺线描述足垫近端月壤的滑移线，采用直线描述远端月壤的滑移线．根据极限土压力理论计算滑移线上微元体的竖向月壤压力和水平向月壤压力，沿曲线积分求得滑移面上月壤的极限承载力，并考虑足垫底部摩擦力及惯性力等因素，建立了足垫-月壤水平拖曳力学模型．对该过程的离散元仿真结果验证了上述理论模型，理论与仿真结果的误差小于9%．1 足垫-月壤水平拖曳力学模型1.1　足垫形状月基装备足垫沿用嫦娥三号着陆腿足垫的结构形式，其底部采用圆形平底设计，如图1所示．10.13245/j.hust.240683.F001图1月基装备足垫形状探月工程中，月基装备包含着陆器、巡视器、功能舱、特种机器人、能源站等多种特种装备，各种装备的足垫尺寸也有所差异．本研究以底部直径6.0 cm的月基装备足垫为例，探究足垫-月壤水平拖曳力学模型．由于各种尺寸的足垫在拖曳过程中的力学原理是相似的，因此本研究揭示的拖曳力学模型也可以推广到其他尺寸的月基装备足垫中．1.2　滑移面月壤阻力1.2.1　滑移面形状水平拖曳开始前，月基装备足垫被埋置在一定深度的月壤中．拖曳过程中，当拖曳力超过静摩擦力时，足垫开始沿拖曳方向在月壤中滑动，周围月壤受到挤压．当壤体应力超过极限承载力时，月壤出现剪切失效，呈现滑移面．图2为月基装备足垫拖曳过程中月壤滑移面的截面图，图中：曲线AB及直线BC为该截面上月壤的滑移线；τAB和σAB分别为滑移线AB上微元段的切应力和正应力；τBC和σBC分别为滑移线BC上微元段的切应力和正应力．10.13245/j.hust.240683.F002图2月壤滑移面在理论模型的建立中，须采用精确的数学表达式描述月壤滑移面的形状和位置．相较于地球土壤，月壤颗粒间的具有更大的滚动摩擦力和滑动摩擦力，在宏观上呈现出更高的内摩擦角．这种高内摩擦角的月壤性质加剧了月壤与足垫间的机械啮合作用，导致近端月壤内部产生压缩应力，提高了该部分月壤的抗剪强度．本研究采用Osman对数螺线[18]描述近端月壤的滑移线AB，即考虑了这种强度提升作用．对于滑移线上远离足垫的月壤而言，这种强度提升作用可忽略不计，因此采用直线描述远端月壤滑移线BC．确定滑移线形状后，分别针对曲线段和直线段计算滑移线上月壤的法向阻力和切向阻力，并最终求得整个滑移面上的月壤阻力．1.2.2　曲线段滑移线阻力曲线段AB采用的对数螺线为[18]ρ=r0ecot(2β)θs，式中：r0和cot(2β)为对数螺线的形状参数，r0为线段OA长度，β=π/4-φ/2为滑移面的滑移角，φ为月壤内摩擦角；ρ和θs为对数螺线上一点的坐标值．对曲线段AB划分微元段，根据月壤极限承载力理论和被动土压力理论[19]，微元段上的切应力和正应力分别表示为τAB=Mr0ecot(2β)θssin(α+β-θ1)+N;      σAB=Mcotφr0ecot(2β)θssin(α+β-θ1)+(N-c)cotφ,式中：c为月壤的黏聚力；α为线段OA和OB的夹角；θ1为微元段与原点的连线与线段OA的夹角；M和N为计算参数，M=14(Kp+1)γsin 2φ,N=12cKpsin 2φ+12c(1+cos 2φ),其中，γ为月壤容重，Kp为月壤压力常数[19]．曲线段AB在水平方向的月壤阻力合力FAB可以通过对τAB和σAB在水平方向上的分量沿曲线AB进行积分计算，即FAB=∫0ατABcosπ4-φ2-α+θ+σABsinπ4-φ2-α+θr0eθstan φdθ. (1)1.2.3　直线段滑移线阻力采用与曲线段相同的计算方法，对直线段BC划分微元段．微元段上的切应力和正应力分别为τBC=Mz+N;σBC=Mzcot φ+(N-c)cot φ,式中z为微元段的深度．直线段BC在水平方向的月壤阻力合力FBC可以通过对τBC和σBC在水平方向上的分量沿直线BC进行积分计算，即FBC=(Mh/2+N)cos βhr1+(Mh/2+N-c)cot φsin βhr1, (2)式中：r1为线段OB的长度；h为B点的深度．1.3　底部摩擦力在水平拖曳过程中，月基装备足垫挤压近端月壤，月壤也会对足垫产生反力作用．这种现象在竖直方向上表现为足垫受到的轴力．由于轴力的存在，足垫底部与月壤产生了滑动摩擦力．钟世英[20]对这种轴力作用进行了理论和试验研究，对于底部半径r=2.65 cm的足垫模型，轴力Fn与拖曳速度之间的关系近似表示为{Fn}N={λv0}N+8.5，式中：λ为轴力传递系数；v0为拖曳速度．考虑到足垫尺寸对轴力的放大作用，足垫在拖曳过程中受到的底部摩擦力可以表示为Ff=Fn(R2/r2)tan φf，(3)式中：R为本研究月基装备足垫的底部半径；tanφf为月壤的残余摩擦角．1.4　惯性力在静力学问题向动力学问题转化的过程中，达朗贝尔原理提供了一种简便而准确的思路，即通过系统的惯性力考虑动力学作用．本研究中，惯性力的来源是受力月壤．将月壤滑移体简化为如图3所示的平行四边形OADC，取平行于线段OA的斜向微元段dx，dx与C点的距离为x．10.13245/j.hust.240683.F003图3惯性力计算示意图当足垫以速度v0在月壤中滑动时，位于O点附近的月壤速度也为v0，而足垫远端(C点附近)的月壤速度为0．将C点到O点的速度变化趋势视为线性递增，则微元段dx的速度为vmic=x2r1cos βv0．因此，微元段dx的惯性力为dFI=x4r1cos βv02dm,式中dm为微元段的质量．对微元段上的惯性力沿CO积分，得到滑移区月壤的惯性力合力为FI=14r12ρmv02sin(2β)sin2β，(4)式中ρm为月壤密度．1.5　水平拖曳力学模型月基装备足垫水平拖曳力(F)的主要组成有式(1)和(2)所示的滑移面月壤阻力、式(3)所示的足垫底部摩擦力和式(4)所示的滑移区月壤惯性力，并以面积相等的原则将滑移区月壤的宽度简化为b=πR．于是，得到完整的月基装备足垫-月壤水平拖曳接触力学模型F=(FAB+FBC+Ff+FI)b．2 拖曳过程离散元建模2.1　足垫物理模型足垫物理模型在三维绘图软件SolidWorks中绘制，其结构形式如图1所示．材料选取月基装备着陆缓冲器常用的超硬铝合金7A09．该材料的特点是强度高、无磁性、耐热性好，且具有良好的塑性．由于其优秀的综合性能，常被用于制造承载较大的月基装备零件．该材料主要物理和机械性能指标如下[21]：强度≥530 MPa，密度为2.80 g/cm3，断裂伸长率≥6%，弹性模量为71 GPa，剪切模量为26.5 GPa．2.2　模拟月壤离散元模型离散元方法最早由Cundall和Strack[22]提出，并将其应用在散粒力学行为的模拟中．离散元方法分析颗粒运动问题的基本思路是将研究对象视为具有特定力学性质元素的集合，对每个元素采用中心差分方法求解运动方程，并计算相邻颗粒的运动及受力关系，最终得到研究对象整体的运动形态．本研究采用离散元软件EDEM对足垫的水平拖曳过程进行仿真分析．月壤是由不同形状、不同粒径的月壤颗粒组成的．根据嫦娥五号取回的CE5C0400号月壤样本[23]的激光衍射图像分析结果，在离散元软件中采用单球模型、双球模型和四球模型模拟典型的月壤颗粒形状，如图4所示．10.13245/j.hust.240683.F004图4月壤颗粒形状及离散元模型考虑到计算机算力的限制，离散元建模时往往不能完全模拟微小的颗粒．但当颗粒过大时，会产生巨大粒径效应，导致仿真的精确度下降．Fuglsang和Ovesen[24]指出，当圆形机具的尺寸大于平均粒径的35倍时，巨大粒径效应的不利影响可以忽略．基于以上考量，并参考月壤及月壤模拟物粒径分布的趋势，仿真中采用的模拟月壤颗粒模型粒径分布为：0.6~0.8 mm粒径的月壤颗粒占10%，0.8~1.0 mm占20%，1.0~1.1 mm占45%，1.1~1.2 mm占25%．本仿真中，模拟月壤颗粒的平均粒径为1.01 mm，足垫底部尺寸为模拟月壤平均粒径的39.6倍，因此巨大粒径效应可以忽略不计．此外，模拟月壤及模拟月壤颗粒的物理参数如下：内摩擦角为43.6°，黏聚力为1.6 kPa，泊松比为0.33，弹性模量为10.1 MPa．2.3　拖曳过程仿真在离散元仿真软件EDEM中将足垫埋置在0.85 cm深度的模拟月壤中，并使其按照设定速度(0.1，0.2，0.3，0.4和0.5 m/s)水平移动．拖曳过程中，月基装备足垫与模拟月壤的接触如图5所示．10.13245/j.hust.240683.F005图5拖曳过程仿真示意图3 仿真结果验证与分析3.1　月壤滑移面形状分析模拟月壤颗粒的移动速度是观察模拟月壤颗粒运动情况及受力模拟月壤整体滑移情况最直观的指标．图6为足垫在模拟月壤中以0.1 m/s的速度进行水平拖曳时，纵向中心截面上模拟月壤颗粒的移动速度分布情况．可以看出：超过最大移动速度20%以上的模拟月壤颗粒大致构成了模型假设中由对数螺线和直线所描述的模拟月壤滑移线，如图6中曲线AB和直线BC所示．因此，足垫-模拟月壤水平拖曳的仿真结果验证了月壤滑移面假定的合理性．10.13245/j.hust.240683.F006图6模拟月壤颗粒移动速度分布情况(色标单位：m/s)3.2　月壤水平阻力分析以0.1 m/s的水平拖曳工况为例，对通过理论计算得到的月壤阻力与通过离散元仿真得到的模拟月壤阻力进行对比，结果如图7所示．10.13245/j.hust.240683.F007图7月壤水平阻力理论与仿真对比由图7可以看出：仿真中，足垫受到的模拟月壤水平阻力在开始阶段快速上升，随后达到并保持相对稳定的运动状态．当拖曳速度同为0.1 m/s时，仿真中，稳定运动阶段模拟月壤的平均阻力为3.51 N；理论模型中，月壤的阻力为3.95 N，理论模型与仿真的误差为12.5%．理论计算中，主要参数的取值为：R=30 mm，r0=10 mm，h=6.4 mm，ρm=1 320 g/cm3，φ=43°，c=1.6 kPa．3.3　结果分析针对速度为0.2，0.3，0.4和0.5 m/s的水平拖曳过程，同样开展了离散元仿真分析和理论对比研究．各个工况的模拟月壤水平阻力仿真值和理论值的对比情况如表1所示，仿真时间步均为1.80×10-5 s．10.13245/j.hust.240683.T001表1仿真与理论结果对比速度/(m∙s-1)仿真阻力/N理论阻力/N误差/%仿真时间/s0.13.513.9512.540.23.794.2311.620.34.304.566.0510.44.974.881.8110.55.895.2111.51根据表1的对比结果进行计算可知，所有工况下的平均误差为8.70%．本研究提出的月基装备足垫-月壤水平拖曳理论模型和仿真结果拟合得较好．由于本研究是探究月基装备行走过程中足垫与月壤的水平滑移，0.5 m/s以下的拖曳速度已经可以覆盖所有的月基装备行走工况．仿真中模拟月壤水平阻力的曲线呈一条斜率略微增加的曲线，而本研究提出的理论曲线更接近一条斜率恒定的直线．这可能说明，拖曳速度对月壤阻力的影响是非线性的，但本研究提出的理论模型对速度影响的考虑还不够全面．除了本研究考虑的足垫底部摩擦力和惯性力受到拖曳速度的影响，还有一些其他因素受速度的影响较大．例如，由于拖曳速度增加，月壤受到的竖向挤压增大，月壤滑移面上的摩擦力也会增大．这些因素可以在将来更进一步的研究中予以考虑．4 结论本研究对月基装备足垫-月壤的水平拖曳过程中的动力学问题展开了理论和仿真研究．基于月壤极限承载力理论，建立了月基装备足垫-月壤水平拖曳力学模型，用以计算各种拖曳速度下的月壤水平阻力．开展了基于离散元方法的模拟月壤水平拖曳仿真研究，并将仿真结果与理论计算结果进行了对比，研究结果如下．a．本研究所提出的月基装备足垫-月壤水平拖曳力学模型具有较高的准确性，月壤水平阻力的仿真结果与理论计算结果符合较好，平均误差仅为8.7%；b．月壤滑移面形状的假定在仿真中得到了验证．这样的滑移面假定同样适用于月壤推铲、挖掘等月壤作业过程；c．研究成果可以与足垫-月壤竖向刺入力学的研究相结合，探究月基装备行走过程中足垫在月壤上的滑移、沉陷等过程，综合研究月基装备的行走稳定性和安全性．