模拟月壤微观结构形态的定量化研究新方法

蒋明镜，奚邦禄，李立青，戴永生

（1.同济大学 土木工程防灾国家重点实验室，上海200092；2.同济大学 岩土及地下工程教育部重点实验室，上海200092）

已有资料表明，砂土的微观结构形态会对其宏观力学响应（如剪切波速、抗剪强度等）产生显著影响，现已引起众多岩土工作者的关注［1-3］.目前，砂土微观结构参数的统计分析方法随测试设备和技术的更新——如电镜扫描［4-6］（scanning electron microscope，SEM）、X射线［7］、PIV（particle image velocimetry）［8］技术、CT技术［9］以及图像分析系统［10］等的运用——已取得很大进步，这对砂土宏观力学性能机理的研究起到很大推动作用.然而，砂土微观结构形态的定量化研究仍以电镜扫描获得高倍清晰照片，再利用图像分析系统或图像处理软件人工分析为主.对于颗粒形态明显、未胶结成块且粒径差别较大的砂土，若能采取简易方法测取所需的微观结构参数，对砂土宏观力学机理研究的发展意义重大.从土力学的角度来看，相对于大量颗粒集合体构成的土体而言，“微观”尺度是指单颗粒的尺度，因此“微观”所对应的绝对尺寸随颗粒大小而变化.

月壤根据级配曲线可归为粉质砂土［11］.本文作者以火山灰为原材料，在文献［12-13］提供的部分月壤样品级配曲线范围内设定几种不同的级配曲线，并对不同样本进行相应的力学特性测试.以真实月壤强度指标为目标，经过10余次的级配调整，研制出同济一号（简称TJ-1）模拟月壤［14-15］.图1为阿波罗11号登月点月壤样品和TJ-1模拟月壤的颗粒级配［14-15］.图中D表示颗粒直径.由图可知，TJ-1模拟月壤在月壤的颗粒级配范围内，其颗粒形态明显，粒径差别较大，符合本文微观参数研究的需求.

图1 TJ-1及部分真实月壤级配曲线［14-15］Fig.1 Gradations of TJ-1and lunar regolith［14-15］

首先通过电镜扫描获取TJ-1模拟月壤各粒组的高倍清晰照片，而后利用AutoCAD和MATLAB软件获取各粒组颗粒的微观结构参数.考虑到TJ-1模拟月壤由7种粒径的颗粒组成，且不同粒组微观结构形态也不尽相同，为此首次尝试采用离散元数值试验近似确定砂样中各粒组所受接触力占试样总接触力的权重，并与各粒组形态参数的均值相结合，采用加权平均的方法获取TJ-1模拟月壤形态参数的代表值.该方法简单易行，不失为一种确定含多种粒径砂土微观结构形态参数的新尝试.

1 微观结构参数及测试方法

1.1 微观结构参数简介

目前，描述土体微观结构形态的参数较多，在此仅介绍部分参数及表达式.表1列出了某典型颗粒基本尺寸参数及描述.图2给出了相应的示意图.

表1 颗粒基本尺寸参数Tab.1 Basic parameters of a particle

颗粒的扁平度即长短轴比例Am为

形状系数Fi为

式中：C为与颗粒等面积圆的周长，即C=πDp.

组构系数［7］FF为

图2 某典型砂土颗粒的基本尺寸参数Fig.2 Basic parameters of a particle

颗粒的粗糙度CR由颗粒等效球（或圆）和颗粒实际封闭曲面（边界轮廓线）间半径均方差来度量［16］.粗糙度的定义按字面上说，对于三维物体即是表面粗糙的程度，对于二维图形即是边界线的曲折程度.本文定义的粗糙度值和圆度定义是一样的，其物理意义是把粗糙度看成是颗粒投影轮廓的实际周长与等效圆的周长之差.对于图2所示的颗粒，其粗糙度CR可由下式表示：

需要指出的是，式（1）～（4）是针对单个颗粒的计算分析，若某粒组中颗粒的形态差别不大，则可选取一定数目的颗粒测试相应的参量，而后取其平均值作为该粒组的形态参量值.但对于由不同粒组组成且级配跨度大的砂样（如TJ-1模拟月壤），微观结构参数代表值的测试方法目前尚未见报道，为此本文拟采用离散元数值试验与单一粒组颗粒形态参数测试相结合的方法进行初步探讨.

1.2 测试方法及过程

刘清秉等［1］指出采用25个颗粒与100个颗粒测试同一粒组砂土的微观结构参数时结果偏差不大；刘瑜等［2］测试砂土的粗糙度参数时选用颗粒数目的平均值为18.上述研究结果表明，在进行单一粒组颗粒微观结构参数的测试时随机选取20个颗粒组成代表性样品基本可以满足要求.本文针对组成TJ-1模拟月壤的7种粒组，采用如下思路测试各粒组的形态参数：首先通过筛分试验获得各粒组的颗粒；然后采用电镜扫描与AutoCAD基本功能相结合的方法，测试同一粒组中典型颗粒的长轴、短轴、周长和面积；而后采用式（1）～（3）获取单个颗粒的长短轴比例、形状系数及组构系数等，取其平均值作为该粒组的微观结构参数（本文仅给出长短轴比例和粗糙度的分析结果）.具体实施步骤如下：

（1）将TJ-1模拟月壤筛分成粒径范围分别为≤0.075mm，＞0.075～0.100mm，＞0.10～0.25 mm，＞0.25～0.50mm，＞0.5～1.0mm，＞1～2 mm，＞2～5mm等7种粒组，清水冲洗、烘干后选取代表性颗粒，对于＞2～5mm粒组的颗粒采用高倍数码相机拍摄，其他粒组采用型号为S-2360N的扫描电镜获取高倍清晰照片.为降低人为主观性造成的试验误差，文中每种粒组选取20个左右的代表性颗粒进行试验.

（2）将扫描和拍摄得到的照片导入AutoCAD，局部放大，保证颗粒外轮廓清晰可见.然后在用多段线描绘颗粒轮廓时，采用尽可能短的线与之贴合.绘制完成后发现，多段线描绘的外轮廓与真实轮廓基本重合，完全能够满足精度要求.

（3）先根据颗粒的形状粗略估计其长轴的大小，据此绘制一个直径相当的圆，然后调整该圆的大小，使其刚好能将颗粒整个包围起来.如图2所示的最小外接圆，其直径即为颗粒的长轴.而后在长轴的中点画一条垂直于长轴并与轮廓线相交的线段，作为该颗粒的短轴.

（4）利用AutoCAD工具栏中的工具→查询→面积/距离命令输出并记录颗粒的周长、面积、长轴和短轴；采用式（1）～（3）获取相应的微观结构参数以及颗粒的等效颗粒直径.

（5）重复步骤（2）～（4）获取该粒组不同颗粒的微观结构参数，取其平均值得到该粒组相应的微观结构参量.

粗糙度CR是在步骤（2）的基础上，通过软件MATLAB中的图像处理功能及自编的MATLAB程序计算获得，具体过程如下：将步骤（2）所得的颗粒轮廓线填充后储存为图片格式（.bmp）；采用MATLAB软件的图像处理功能将图片转化为黑白二元图，利用Canny算子得到该颗粒的轮廓曲线；而后采用自编的MATLAB程序按逆时针方向获取轮廓线上各点坐标并转化为坐标原点在颗粒形心的极坐标数列Zi（Ri，θi），求得极坐标幅值的平均值R，采用式（4）即可求得该颗粒的粗糙度.以圆形、正方形和长短轴比例Am=1.5，2.0，3.0的椭圆进行粗糙度计算精度的验证，具体流程见图3（以圆形、正方形为例），计算精度见表2.由表2可知，采用自编MATLAB程序测得圆形的粗糙度与理论值吻合良好，误差不超过0.000 1%，正方形的粗糙度与理论值的误差也仅为3.31%，同时可以发现随长短轴比例（扁平度）的增加，椭圆形颗粒的粗糙度呈增大趋势，故而式（4）在一定程度上包含有长短轴比例对粗糙度的影响.因此采用本方法获取颗粒表面粗糙度CR是可行的.

表2 颗粒粗糙度的计算误差Tab.2 Computational errors of particle roughness

图3 颗粒粗糙度的获取过程Fig.3 Procedure of getting the roughness of particles

2 试验结果及分析

2.1 TJ-1模拟月壤的颗粒形态及测试结果分析

真实月壤颗粒形态较为多变，从圆球状、椭球状到极端棱角状都有，但以长条状、次棱角状和棱角状较为常见［17］.图4为TJ-1模拟月壤部分粒组及典型颗粒的电镜扫描图像.由图4a～d可知，TJ-1模拟月壤同一粒径组中，颗粒形状差别很大，从长条状到接近圆形都有分布.颗粒多数呈长条状，存在棱角、钩角、锯齿等不规则结构，与月壤的颗粒形态较为相似［17-18］，且长短轴比例Am一般不大于4.这样的颗粒形态致使TJ-1模拟月壤遭受外力作用时能够比较容易地相互啮合形成受力骨架，抵抗一定的外力和变形.当颗粒粒径较大时（图4d），部分颗粒带有蜂窝状气孔结构，而当颗粒粒径较小时（图4a～c），则没有蜂窝状气孔结构，仅由细小的火山灰颗粒组成.从典型颗粒的电镜扫描图像（图4e）可知，颗粒表面凹凸不平、较为粗糙.

将SEM扫描和高清相机拍照获得的图像导入到AutoCAD中进行描绘处理后获取各粒组颗粒的有效图像数量如表3所示.表4和表5分别为采用1.2节所述方法分析计算得到的TJ-1模拟月壤各粒组的长短轴比例及粗糙度的试验结果.图5和图6为两种形态参数的分布情况.

图4 模拟月壤部分粒组及典型颗粒的电镜扫描图像Fig.4 SEM images of TJ-1lunar soil simulant

表3 TJ-1模拟月壤各粒组扫描获得的有效颗粒数目Tab.3 Particle number of each fraction in SEM tests

图5 TJ-1模拟月壤各粒组颗粒的长短轴比例分布Fig.5 Aspect ratio of each fraction of TJ-1lunar soil simulant

由表4、表5及图5、图6的计算分析结果可知：TJ-1模拟月壤颗粒的长短轴比例Am在1.149～3.283之间，而各粒组的平均值变化范围较小，仅在1.615～1.815之间波动；粗糙度参量CR在0.025 5～0.156 7之间，各粒组的平均值变化范围为0.068～0.082.

2.2 TJ-1模拟月壤形状参数的确定

为尽量减小粒径分布对砂样宏观力学响应的影响，文献［1-2］均采用单一粒组的砂样研究微观结构形态对其力学性质（如剪切波速、抗剪强度等）的影响.然而自然界中砂土是由不同粒组的颗粒组成，如何测试多粒组组成的砂土微观结构参数，进而确定微观结构对宏观力学性能影响的研究尚未见报道，这也是土体微观结构研究方面的难题.考虑到砂土受力过程中同一粒组颗粒所受接触力差别较小［3］，不同粒组的接触力大小与粒径直接相关.故而可考虑采用各粒组颗粒接触力之和占试样中颗粒间总接触力的权重作为相应各粒组微观结构参数的权重.然而，物理试验中颗粒间的接触力无法测得，为此，本文尝试采用离散单元法进行TJ-1

表4 TJ-1模拟月壤长短轴比例的统计结果Tab.4 Calculated aspect ratio of TJ-1lunar soil simulant

表5 TJ-1模拟月壤粗糙度参数的统计结果Tab.5 Calculated roughness of TJ-1lunar soil simulant

图6 TJ-1模拟月壤各粒组颗粒的粗糙度分布Fig.6 Roughness of each fraction of TJ-1lunar soil simulant

模拟月壤离散元试样的等向固结数值试验，计算试样中各粒组颗粒的接触力总和，求得各粒组的接触力占试样总接触力的权重，并将该粒组接触力的权重作为相应的微观结构参数的权重，采用加权平均求得该结构参数的代表值.该方法为测定单一粒组微观结构参数的常规测试技术与离散元数值分析相结合的新型测试方法，是测定由不同粒组组成砂样的微观形态参数的新尝试.具体实施过程如下：

（1）首先采用1.2节中的方法测得TJ-1模拟月壤各粒组的微观结构参数.

（2）采用分层欠压法［19］生成TJ-1模拟月壤的离散元试样，其级配曲线如图7所示.图中各粒径为1.2节测得的各粒组等效颗粒直径平均值，各粒径所对应的质量百分比见文献［12］.离散元试样的颗粒总数为24 000，高宽比为2：1.通过二维离散元单元试验模拟真实月壤力学特性，选取平面目标孔隙比e0=0.22.若与实际模拟月壤孔隙比相对应，可以通过以下公式［20］类推：

式中：b，h分别为试样尺寸的宽度与高度；N为总颗粒数目；ri为颗粒半径；r-为平均粒径.

图7 确定接触力权重时TJ-1模拟月壤的颗粒级配Fig.7 Distribution of grain size of TJ-1lunar soil simulant in DEM analyses

在离散元模拟过程中以圆形颗粒为主，难以模拟实际颗粒的形态.因为考虑实际颗粒形态，需要用复杂的函数取逼近颗粒形状，效率非常低，而且在加载过程中，颗粒形状也会改变.本文研究不同粒径颗粒在加载过程中承受荷载的比值，因此主要针对颗粒大小，即颗粒级配这一主要特征.颗粒形态虽然与真实图片有差距，但还是能够满足本文研究的需要.

离散元试样首先在12.5kPa下单向固结，而后在100kPa下等向固结稳定（图8），随后计算输出试样中各粒组颗粒所受接触力的总和以及试样中所有颗粒接触力的总和，进而求得每种粒径的接触力占总接触力的权重.分层欠压法生成离散元试样的基本原理以及试样生成时每种粒径所对应颗粒数目的确定方法已在文献［19］中详细叙述，在此直接给出数值试验结果，见表6.

图8 TJ-1模拟月壤离散元试样等向固结示意图Fig.8 Schematic diagram of isotropic consolidation of DEM sample of TJ-1lunar soil simulant

表6 TJ-1模拟月壤各粒组所受力的权重Tab.6 Weight of contact force of each fraction in sample

（3）将表6中的计算结果与各粒组的微观结构参数（见表4和表5）相结合，采用加权平均，即可求得TJ-1模拟月壤相应的微观结构参数.

由表6可知，离散元试样等向固结完成后，各粒组颗粒所受接触力的权重差别明显，粒径越小的粒组在离散元试样中数目越多，相应接触力的权重越大，从而能够在受力骨架中占主导作用；而＞1～2 mm和＞2～5mm粒组的颗粒在离散元试样中数目很小，相应的接触力权重也很小，仅为0.483%～1.558%.综合表4～6的试验结果，采用加权平均求得TJ-1模拟月壤长短轴比例和粗糙度参数分别为1.711和0.077.上述微观结构参数可直接与本团队改进后的NS2D离散元程序相结合用于TJ-1模拟月壤力学性能的研究，采用该程序探索砂土的非共轴微观机制也是本团队研究内容之一.

3 结语

传统的微观结构测试主要针对单一粒组进行，本文针对TJ-1模拟月壤这种多粒径的颗粒群，提出一种室内试验与离散元数值模拟相结合的方法用以测定由不同粒组组成的砂样微观结构参数.首先通过电镜扫描获取TJ-1模拟月壤各粒组的高清照片，利用AutoCAD和MATLAB软件方便快捷地获取各粒组颗粒的长短轴比例及粗糙度等微观结构参量；而后针对目前仅能测定单一粒组砂样的微观结构参数的缺陷，通过离散元数值试验近似确定TJ-1模拟月壤中各粒组颗粒所受接触力占所有颗粒接触力的权重，并与各粒组的形态参数相结合采用加权平均来确定其形态参数的代表值.

不同土体有各自的代表性形态参数，该参数由颗粒形态和颗粒级配共同控制.采用本文提出的加权方法获得的形态参数能考虑颗粒形态和颗粒级配的综合影响，具有更好的代表性，这也极大地方便与其他土体的对比，以确定影响TJ-1力学特性的关键因素.新方法简单易行，所得参数可直接应用到已有离散元程序中进行土体力学性质的研究，不失为确定砂土微观结构形态参数的新尝试.

［1］ 刘青秉，项伟，Budhu M，等.砂土颗粒形状量化及其对力学指标的影响分析［J］.岩土力学，2011，32（增1）：190.LIU Qingbing，XIANG Wei，Budhu M，et al.Study of particle shape quantification and effect on mechanical property of sand［J］.Rock and Soil Mechanics，2011，32（Suppl.1）：190.

［2］ 刘瑜，夏唐代.砂土颗粒粗糙度对剪切波速影响的试验研究［J］.岩土工程学报，2011，33（2）：285.LIU Yu，XIA Tangdai.Experimental study on influence of particle roughness on shear wave velocity of sand ［J］.Chinese Journal of Geotechnical Engineering，2011，33（2）：285.

［3］ JIANG Mingjing，Leroueil S，ZHU Hehua，et al.Twodimensional discrete element theory for rough particles ［J］.International Journal of Geomechanics，2009，9（1）：20.

［4］ 唐朝生，施斌，王宝军.基于SEM土体微观结构研究中的影响因素分析［J］.岩土工程学报，2008，30（4）：560.TANG Chaosheng，SHI Bin，WANG Baojun.Factors affecting analysis of soil microstructure using SEM ［J］.Chinese Journal of Geotechnical Engineering，2008，30（4）：560.

［5］ 王慧妮，倪万魁.基于计算机X射线断层术与扫描电镜图像的黄土微结构定量分析［J］.岩土力学，2012，33（1）：243.WANG Huini，NI Wankui.Quantitative analysis of loess microstructure based on CT and SEM images［J］.Rock and Soil Mechanics，2012，33（1）：243.

［6］ 李丽华，唐辉明，刘数华.月壤及模拟月壤微观结构的研究［J］.岩土力学，2012，33（1）：31.LI Lihua，TANG Huiming，LIU Shuhua.Microstructure of lunar soil and lunar soil simulant［J］.Rock and Soil Mechanics，2012，33（1）：31.

［7］ Alshibli K A，Alsaleh M I.Characterizing surface roughness and shape of sands using digital microscopy ［J］.Journal of Computing in Civil Engineering，2004，18（1）：36.

［8］ White D J，Take W A，Bolton M D.Soil deformation measurement using particle image velocimetry （PIV）and photogrammetry［J］.Géotechnique，2003，53（7）：619.

［9］ 孙红，葛修润，牛富俊，等.上海粉质粘土的三轴CT实时细观试验［J］.岩土力学与工程学报，2005，24（24）：4559.SUN Hong，GE Xiurun，NIU Fujun，et al.Real-time CT mesotesting on Shanghai silty clay subjected to triaxial loading ［J］.Chinese Journal of Rock Mechanics and Engineering，2005，24（24）：4559.

［10］ SHI Bin，LI Shenglin，Tolkachev M.Quantitative approach on SEM images of microstructure of clay soils ［J］.Science in China：Series B，1995，38（6）：741.

［11］ Heiken G，Vaniman D，French B M.Lunar sourcebook：A user’s guide to the moon ［M］.New York：Cambridge University Press，1991.

［12］ Chang C S，Hicher P Y.Model for granular materials with surface energy forces ［J］.Journal of Aerospace Engineering，2009，22（1）：43.

［13］ Mitchell J K，Bromwell L G，Carrier W D，et al.Soil mechanical properties at the Apollo 14 site［J］.Journal of Geophysical Research，1972，77（29）：5641.

［14］ 蒋明镜，李立青.TJ-1模拟月壤的研制［J］.岩土工程学报，2011，33（2）：209.JIANG Mingjing，LI Liqing.TJ-1 lunar soil simulant ［J］.Chinese Journal of Geotechnical Engineering，2011，33（2）：209.

［15］ JIANG Mingjing，LI Liqing，SUN Yugang.Properties of TJ-1 lunar soil simulant［J］.Journal of Aerospace Engineering，2012，25（3）：463.

［16］ Carrier W D，Olhoeft G R，Mendell W.Physical properties of the lunar surface，in lunarbook［M］.Cambridge：Cambridge University Press，1991.

［17］ 欧阳自远.月球科学概论［M］.北京：中国宇航出版社，2005.OUYANG Ziyuan.Introduction to lunar science［M］.Beijing：China Astronautic Publishing House，2005.

［18］ 胡荣泽，刘森英.颗粒的表面粗糙度［J］.粉体技术，1998，4（3）：1.HU Rongze，LIU Senying.Surface roundness of particles［J］.Powder Science and Technology，1998，4（3）：1.

［19］ JIANG Mingjing，Leroueil S，Konrad J M.An efficient technique for generating homogeneous specimens for DEM studies［J］.Computers and Geotechanics，2003，30（5）：579.

［20］ 周健，张刚，孔戈.渗流的颗粒流细观模拟［J］.水利学报，2006，37（1）：28.ZHOU Jian，ZHANG Gang，KONG Ge.Meso-mechanics simulation of seepage with particle flow code［J］.Journal of Hydraulic Engineering，2006，37（1）：28.