随着嫦娥五号成功从月表取样返回，圆满完成“绕、落、回”三步走战略，我国的探月工程已进入以月面资源原位利用和建设科研站基本型为主要目标的新阶段．即将实施的探月工程四期任务将开展大规模的月球综合地质勘查、月壤钻探及取样和月面资源原位利用等工程任务．月壤是月球探测的主要研究对象和载体，系统研究并掌握月壤的基本物理力学性质是探月工程计划的首要任务，也是成功实施探月工程任务的重要保障．月壤是在陨石撞击、太阳风粒子轰击及月面环境影响的共同作用下形成的，其独特的形成过程导致月壤具有特殊的颗粒形态，从而使其物理力学特性与地球土壤显著不同．通过月面原位勘探及阿波罗任务带回的月壤样品地面试验研究加深了人们对于月壤力学性质的认识．另外，由于月壤样品的稀缺性，国内外学者大多都采用与月壤具有相似物理力学特性的模拟月壤进行相关实验研究．出于不同的研究目的，各国学者已经研制了数十种模拟月壤[1]．包括美国的MLS-1，JSC-1/A和GRC-1，日本的FJS-1及我国的CAS-1和TJ-1等模拟月壤．通过对这些模拟月壤开展室内试验研究，极大地丰富了人们对于月壤力学特性及月壤与结构物相互作用机制的认识．除了在地面进行室内试验研究外，采用离散单元法从颗粒属性和颗粒间相互作用效应出发，通过数值模拟研究月壤的力学特性同样是一种有效可行的方法．早期的离散元数值模拟通常将月壤颗粒简化为理想的圆或球颗粒，忽略了其复杂形状的影响．Hasan等[2]通过PFC3D研究了JSC-1A模拟月壤在不同围压和密实度下的峰值和临界摩擦角的变化特点．Li等[3]研究了QH-E模拟月壤在恒定平均主应力路径下的应变软化和剪胀特性并探讨了细观参数的影响．然而，已有的实验和数值模拟研究[4-6]均表明，颗粒材料的形状特征会显著改变其力学特性．月壤具有特殊的颗粒形态[7]，形状差异大，从球形到尖锐棱角状等不规则形状都有分布，且表面凹凸不平、粗糙度较大．因此，Jiang等[8]通过在接触模型中添加滚动阻力，部分考虑了颗粒的形状效应．Katagiri等[9]通过CT扫描获取了74个月壤颗粒的三维形状，并结合离散元法研究了形状对休止角和简单剪切行为的影响．Khademian等[10]研究了四种特定形状对月壤颗粒在不同重力条件下休止角大小的影响．目前，针对月壤及模拟月壤的力学特性，许多学者已经进行了大量研究，然而，基于月壤及模拟月壤颗粒实际形状的离散元数值模拟研究还相对较少，颗粒形状对其宏细观力学特性的影响研究尚不系统和深入．因此，本研究以IRSM-1模拟月壤[11]颗粒为研究对象，利用图像处理技术对模拟月壤颗粒的形态特征进行定量描述；在此基础上，将形状数据导入PFC2D离散元软件中，实现真实颗粒形状的离散元模拟；建立具有不同整体规则度的模拟月壤数值试样，开展双轴压缩数值试验，系统研究模拟月壤颗粒形状对宏细观力学行为的影响机制．1 模拟月壤颗粒形态特征1.1　IRSM-1模拟月壤颗粒轮廓提取IRSM-1模拟月壤[11]是以玄武质火山灰和钛铁矿为原料，经过高温烧结和冲击破碎等工艺制备而成的，属于高钛模拟月壤．其相对密度为3.26，黏聚力和内摩擦角分别为0.90 kPa和44.34°．IRSM-1模拟月壤较好地还原了真实月壤的物理力学特性及复杂的颗粒形态特征，可以用于月壤性质的研究．鉴于该模拟月壤颗粒较小，本研究通过显微镜和图像处理技术获取了IRSM-1模拟月壤颗粒的轮廓数据，流程如图1所示．具体的操作步骤如下：首先采用工业显微镜拍摄模拟月壤颗粒的高清数字图像；然后将其转化为灰度图像，并通过中值滤波法进行图像增强；最后，通过对灰度图像进行二值化处理，即可提取颗粒轮廓点的坐标．10.13245/j.hust.240849.F001图1IRSM-1模拟月壤颗粒轮廓点坐标提取流程1.2　模拟月壤颗粒形态特征在获得模拟月壤颗粒的轮廓坐标后，可计算相应的颗粒形态评价指标．选取颗粒的纵横比(δAR)、圆形度(S)和凹凸度(C)作为颗粒的形状参数，如图2所示．具体定义如下：δAR为颗粒的短轴长度Lmin与长轴长度Lmax之比；S为与颗粒面积相同的等效圆的周长2πre与颗粒周长P之比，其中，re为等效圆的半径；C为颗粒面积A与其凸包面积Ac之比．10.13245/j.hust.240849.F002图2颗粒形状参数定义及计算方法示意图为了定量评价IRSM-1模拟月壤颗粒的形态特征，共获取1 250个模拟月壤颗粒的轮廓坐标并计算了相应的形状参数．图3给出了纵横比、圆形度和凹凸度的累计平均值及累积分布曲线图．从图3(a)可以看出：随着统计颗粒数目的增加，模拟月壤颗粒形状参数的平均值逐渐趋于稳定；当统计的颗粒数目达到900时，纵横比、圆形度和凹凸度的累计均值不再变化，因此选取的1 250个颗粒足以反映IRSM-1模拟月壤颗粒的形态特征．参考已有的研究[4-6]，选择累积分布为50%所对应的形状参数值作为模拟月壤颗粒的形态特征值，即纵横比(δAR50)、圆形度(S50)和凹凸度(C50)的特征值分别为0.743，0.755和0.902．10.13245/j.hust.240849.F003图3模拟月壤颗粒的形状参数统计结果由于前述三个形状参数仅反映了模拟月壤颗粒某一方面的形态特征，本研究选用整体规则度(R)作为颗粒形态特征的综合表征[4]，其定义为颗粒的纵横比、圆形度及凹凸度之和的平均值，即R=(δAR+S+C)/3．(1)整体规则度越高(接近1)，表明颗粒越规则(接近圆形)；反之，颗粒越不规则．本研究模拟月壤颗粒的整体规则度RC=0.8，而圆形颗粒的整体规则度RR=1.0．根据Yang等[4]提出的组合形状参数的概念，定义试样的整体规则度取决于对应颗粒材料的整体规则度及其百分比．例如，由50%的模拟月壤颗粒和50%圆形颗粒组成的数值模拟试样C50R50的组合整体规则度为RC50R50=50%RC+50%RR．通过将一定比例(0%，25%，50%，75%，100%)的圆球颗粒随机地替换成具有模拟月壤实际形状的颗粒簇，可以制备得到不同整体规则度的离散元数值试样．表1总结了本研究使用的五种试样的形状参数，表中C0R100和C100R0分别代表全部由圆颗粒和模拟月壤颗粒组成的试样．10.13245/j.hust.240849.T001表1不同模拟月壤颗粒含量试样的形状参数试样纵横比圆形度凹凸度整体规则度C0R1001.0001.0001.0001.00C25R750.9360.9390.9760.95C50R500.8720.8780.9510.90C75R250.8070.8160.9270.85C100R00.7430.7550.9020.802 离散元数值模拟2.1　单一颗粒实际形状模拟采用二维离散元软件PFC2D来分析模拟月壤在双轴压缩条件下的宏细观力学行为．在提取了模拟月壤颗粒轮廓后，将其导入PFC软件中，并采用泡沫填充算法[12]生成具有真实形状的颗粒簇Clump．该算法有ρ和φ两个控制参数，其中：ρ为构成颗粒簇Clump的圆颗粒的最小和最大半径之比；φ为相邻两个泡沫圆颗粒的交角．当φ=0时，两个圆颗粒外部相切；当φ=180°时，一个小圆内切于另一个小圆．不同控制参数组合情况下真实颗粒形状的模拟效果如图4所示．表2给出了不同控制参数组合下颗粒簇面积与颗粒实际面积之比．从图4中可以看出：随着ρ的减小和φ的增加，模拟同一个实际形状所使用的圆颗粒数量增加，模拟的精度也提高．为了综合考虑模拟精度和计算效率，选取ρ=0.1，φ=130°来生成IRSM-1模拟月壤颗粒的实际形状，对应的颗粒簇由37个圆颗粒组成，且面积与颗粒实际面积之比为96%，较好地还原了该颗粒的实际形状．通过该方法一共生成了1 250个不同形状的模拟月壤颗粒簇模板．10.13245/j.hust.240849.F004图4不同控制参数下颗粒簇模拟真实形状效果对比10.13245/j.hust.240849.T002表2不同控制参数组合下颗粒簇面积与颗粒实际面积之比φρ0.050.100.200.3090°94.4190.7686.3475.71110°97.2294.1390.0582.53130°97.8296.0090.9281.66150°98.4797.1293.5289.70%2.2　双轴压缩离散元模型采用以下步骤来生成包含模拟月壤颗粒的双轴压缩试样．首先，在水平和竖直方向生成刚性墙体，并在墙体范围内生成密实均匀的圆形颗粒，级配曲线如图5所示，图中：δu为颗粒的不均匀系数；de为圆颗粒和对应颗粒簇Clump的等效直径，de60，de50和de10分别为累计分布60%，50%和10%所对应的等效直径．然后，按照面积相等的原则，随机挑选一定数量的圆形颗粒，并将其替换为具有实际形状的颗粒簇Clump．颗粒簇的形状是从前述生成的1 250个模板中随机选定的，从而保证了生成的颗粒与真实的模拟月壤颗粒具有相似的形态分布特征．最后，利用伺服机制使试样在给定的围压下进行等向压缩，当模型中的不平衡力与接触力的平均值之比小于1×10-5时，认为试样已达到平衡状态．10.13245/j.hust.240849.F005图5颗粒级配曲线等向压缩时，颗粒间的滑动摩擦系数设置为0，以确保各个不同的试样在轴向加载前具有相同的密实状态．此外，由于颗粒的几何形态特征还包含表面纹理，例如颗粒的粗糙度等，而本研究采用Clump模拟颗粒实际形状的方法还无法实现对颗粒粗糙度的有效表征，因此通过调整颗粒间的摩擦系数来反映表面纹理的影响．生成的离散元模型尺寸约为38 mm×76 mm，共包含3 521个颗粒．不同颗粒簇Clump含量的离散元试样如图6所示．10.13245/j.hust.240849.F006图6不同颗粒簇Clump含量的离散元试样试样固结完成后，通过恒定的应变速率移动上下墙体，同时维持左右墙体的围压不变以模拟双轴剪切试验，直到试样的轴向应变达到15%．为了维持准静态加载条件，设置上下墙体的加载应变率ε˙=0.2 s-1，以确保惯性指数[13]I=ε˙m/p0.001，式中：m为试样颗粒质量；p为围压．由于本研究的重点在于模拟月壤的颗粒形状对其力学特性的影响，因此颗粒间的接触采用简单的线性模型，且颗粒簇Clump和圆颗粒的细观参数保持一致．模型参数如下：颗粒密度为3 100 kg/m3，颗粒间滑动摩擦系数为0.7，颗粒-墙体间滑动摩擦系数为0.1，阻尼系数为0.7，有效模量为1×107 Pa，颗粒法向和切向刚度比为1.0．图7为C100R0试样不同围压下的偏应力-应变曲线，图中：σd为偏应力；ε为轴向应变；σ3为围压．可以看出，在偏应力达到峰值后，试样表现出明显的应变软化特性．10.13245/j.hust.240849.F007图7C100R0试样不同围压下的偏应力-应变曲线此外，根据莫尔-库仑强度准则进行线性拟合后，可计算该试样的黏聚力为1.58 kPa，内摩擦角为43.94°，与室内试验测得的IRSM-1模拟月壤[11]黏聚力(0.9 kPa)和内摩擦角(44.34°)较为符合，也与已有文献中报道的月壤黏聚力(0.8~2.1 kPa)和内摩擦角(37°~45°)[14]较为接近．因此，可以认为本文离散元模型及选取的相关模型参数反映了真实月壤的力学特性，具有一定的合理性．3 离散元模拟结果分析3.1　宏观力学响应3.1.1　应力-应变关系图8展示了在围压为50 kPa的条件下，不同整体规则度试样双轴压缩模拟结果．可以发现，试样的整体规则度对其变形特性、抗剪强度及体变特性都有显著影响．随着整体规则度的下降，即试样中具有模拟月壤实际形状的颗粒数量增多，试样的峰值强度不断提升．例如，全部为实际形状颗粒的试样C100R0的峰值强度约为圆颗粒试样C0R100的两倍．此外，在试验初始加载阶段，曲线近似呈线性增长，曲线的斜率随着试样颗粒不规则度的增加而增大，即试样的整体规则度越低，刚度越大．这种现象可以解释为颗粒形状会影响颗粒之间的细观接触和力链传递．不规则的模拟月壤颗粒间的咬合和互锁作用更为显著，形成的接触点相较于圆形颗粒更多且更稳定，因此从宏观上表现为不规则试样的初始刚度更大．10.13245/j.hust.240849.F008图8不同整体规则度试样双轴压缩模拟结果图8(b)展示了在不同整体规则度下，试样的体应变与轴向应变的关系曲线．从图中可以发现：随着轴向应变的增大，所有试样均表现出先剪缩(体应变为负)后剪胀(体应变为正)的特性；同时，随着试样中不规则颗粒含量的增加，试样的剪胀性逐渐增强．这种现象可能是由于不规则形状的颗粒在绕相邻的颗粒滚动或滑动时会产生更大的相互错动．相比之下，规则的圆颗粒能够以较小的幅度发生滚动和滑动，从而抑制试样的剪胀行为．3.1.2　摩擦角和剪胀角为了定量讨论模拟月壤颗粒形状对其抗剪强度及剪胀特性的影响，计算不同试样的峰值和临界摩擦角，具体为φp=arc sinσ1p-σ3σ1p+σ3; (2)φc=arc sinσ1c-σ3σ1c+σ3, (3)式中：φp和φc分别为峰值和临界摩擦角；σ1p和σ1c分别为峰值和临界轴向应力；σ3为围压，此处取值50 kPa．为避免临界状态阶段应力波动的影响，σ1c取10%~15%轴向应变范围内轴向应力的平均值．图9给出了摩擦角和剪胀角(φd)随试样整体规则度的变化，其中φd=φp-φc．从图中可知：全部由圆颗粒组成的试样C0R100的峰值摩擦角和临界摩擦角分别为33.25°和20.04°，而全部由模拟月壤实际形状颗粒组成的试样C100R0的相应值分别为44.86°和27.88°，考虑模拟月壤真实形状的试样抗剪强度有明显的提升．此外，还可以观察到较小整体规则度的试样具有相对大的剪胀角，这一结论与前述体应变的变化规律相符合．10.13245/j.hust.240849.F009图9摩擦角和剪胀角随试样整体规则度的变化3.2　细观力学特性3.2.1　力学配位数配位数是指试样中某一颗粒与其相邻颗粒的接触数目，即法向接触力大于0的接触数目．它受试样孔隙率、颗粒形状和级配等因素的影响，是定量描述颗粒材料内部结构和稳定性的重要指标．为剔除对颗粒材料力学稳定无贡献的接触，选用平均力学配位数Zm，定义为Zm=2C-N1N-N1-N0，(4)式中：C和N分别为试样中接触和颗粒的总数；N0为周围无激活接触的颗粒数，即悬浮颗粒的数目；N1为周围接触数目为1的颗粒数．平均力学配位数越大，表明试样颗粒之间相互作用更加紧密，颗粒周边孔隙减小；反之，平均力学配位数越低，颗粒周边更容易产生较大的孔隙结构而产生体积膨胀．图10展示了不同整体规则度试样的平均力学配位数与轴向位移关系曲线．可以发现，整体规则度越低，试样的平均力学配位数越大．这是由于模拟月壤颗粒形状复杂，单个颗粒的接触点相较于圆形颗粒更多．具体而言，圆形颗粒试样C0R100的平均力学配位数最小，且在剪切过程中逐渐降低直至趋于稳定值．而包含有模拟月壤颗粒实际形状的试样的平均力学配位数在剪切过程中先增加然后逐渐下降至某一稳定值．这种配位数的变化是由于轴向剪切导致试样内部孔隙被压缩，颗粒间有效接触增多．试样中不规则颗粒的数量越多，配位数增长的幅度越大．然而，随着剪切的继续，应力达到峰值后试样发生剪胀，颗粒间孔隙增多，有效接触减少导致配位数下降．最终，在试样内部颗粒形成稳定的力链结构，平均力学配位数趋于稳定值．10.13245/j.hust.240849.F010图10平均力学配位数与轴向位移关系曲线3.2.2　法向接触力已有研究表明，法向接触力对偏应力张量的影响比切向接触力更为显著[5]，因此本研究主要关注法向接触力在剪切过程中的变化情况．为了反映法向接触力的空间分布各向异性，采用Rothenburg等[15]提出的傅里叶级数拟合方法对其进行描述，具体表达式为f¯n(θ)=f¯0[1+ancos2(θ-θn)]，(5)式中：f¯n(θ)为接触法向落在某角度区间的所有法向接触力的平均值；f¯0为所有接触点的平均法向接触力；an为傅里叶拟合参数，可反映法向接触力各向异性的大小；θ为法向接触力的方向；θn为法向接触力各向异性的主方向[16]．图11给出了不同整体规则度试样在双轴压缩初始、峰值和最终状态(轴向应变为15%)下平均法向接触力(f¯n/f¯0)分布图，图中每一根棱条的长度和指向分别代表平均法向接触力的大小和方向，虚线为式(5)中傅里叶级数的拟合曲线．从图中可以看出：当未发生剪切时，直方图和拟合曲线均近似圆形，表明不同形态特征试样的法向接触力基本表现为相似的各向同性；随着轴向剪切应变的增加，在峰值状态下，直方图和拟合曲线呈“花生状”，平均法向接触力由均匀分布向竖向集中转变，表现出明显的各向异性；当到达最终状态时，法向接触力的分布依旧表现出与峰值状态相同的各向异性特点，但是各向异性程度有所降低．10.13245/j.hust.240849.F011图11不同整体规则度试样在双轴压缩初始、峰值和最终状态下平均接触法向力分布图对比不同试样的法向接触力分布图可以发现：在峰值阶段，整体规则度高的试样法向接触力的各向异性更强，其中圆形颗粒组成的试样C0R100的各向异性明显强于全部由模拟月壤颗粒组成的试样C100R0．这一现象可能的解释为：具有较低整体规则度的C100R0试样中存在更多的有效接触，且这些接触的大小和方向分布较均匀，从而导致峰值阶段的法向接触力的各向异性减弱．这一解释与前述C100R0试样具有最大的平均配位数的结果相符合．整体而言，试样的颗粒形态特征对于法向接触力在剪切过程中的各向异性产生了一定的影响，且这一影响在峰值阶段尤为明显．4 结论本研究基于图像处理技术与离散元方法，将具有IRSM-1模拟月壤实际形状的颗粒簇和圆颗粒按一定的比例混合，制备了五种不同整体规则度的离散元试样，并开展了双轴压缩数值试验．系统地探讨了模拟月壤的颗粒形状对其宏细观力学特性的影响，得出的主要结论如下．a．在宏观力学特性方面，颗粒形状对IRSM-1模拟月壤试样的应力-应变关系和抗剪强度具有重要影响．形态复杂、整体规则度越低的试样，由于颗粒间具有更强的咬合和互锁作用，其初始刚度、抗剪强度、峰值和临界摩擦角均越高，且表现出更显著的应变软化效应．忽略月壤及模拟月壤复杂的颗粒形状、基于理想圆或球颗粒开展的离散元模拟研究可能会低估其工程力学性能．b．在细观力学特性方面，形态复杂的模拟月壤试样具有更大的平均配位数，内部接触更多且更稳定，因此表现出更高的抗剪强度．不规则颗粒导致试样等向压缩后内部存在更多孔隙，在剪切过程中孔隙面积先减少后增加，使平均配位数表现出先增加后减少至趋于稳定的趋势．试样形态越复杂，内部孔隙越多，平均配位数增加的幅度越大．较大的平均配位数导致试样峰值阶段的法向接触力各向异性的减弱．须要注意的是，本研究仅涉及模拟月壤颗粒的二维形状影响，并未充分考虑试样的初始状态、颗粒级配及更真实复杂的三维形态等因素；因此，未来有必要针对上述因素进行更深入的考察，以进一步补充和完善对月壤及模拟月壤力学性质的理解．