基于SEM图像的灰土力学模型研究

岳建伟，陈 颖，孔庆梅，赵丽敏，罗振先

(1. 河南大学土木建筑学院，河南开封 475004； 2. 河南大学开封市不可移动文物修复与安全评价重点实验室，河南开封 475004； 3. 中南大学土木工程学院，湖南长沙 410075 )

0引 言

中国现存的城墙、土遗址和桥梁等不可移动文物建造中均使用了灰土[1-3]，这与灰土结构致密、强度高、与古建筑兼容性好密切相关[4-5]。随着时间的推移，在长期自然和人为因素的影响下，以灰土为主的城墙出现了大量裂隙及表面不同程度的剥落现象，这些裂隙对城墙的稳定性及长期保存产生重要影响。以灰土为主的城墙均为文物，不能随意取土进行力学性能研究[6-7]，相对与现代岩土而言，城墙灰土的力学性能及劣化性能研究相对困难，且受温湿、冻融、盐渍等多种作用影响，进一步增加了城墙灰土研究的难度[8-9]。针对文物保护的特殊要求，探寻城墙灰土力学性能的研究方法是分析灰土劣化和文物本体安全评价关键问题之一。

灰土由石灰、黏土、石英砂等材料组成，各相材料相互接触胶结效应明显，呈现出非均质性和随机性[10-11]。根据灰土的胶结特性可将土体分为土颗粒、胶结物和孔隙三部分[12-13]。胶结物附着在土颗粒接触处的间隙起连接颗粒的作用，其集片状矿物集合体杂乱分布于颗粒之间的空隙之中或呈面-面相接、紧密定向排列于颗粒之间，形成几何形态各异的粒间胶结。因此，从灰土三相介质的力学角度入手，准确把握各相材料组分，尤其是颗粒与胶结物之间的强度差异和相互作用，才能更精确地确定灰土的受力特性以及损伤破坏情况。

土颗粒间的胶结作用是影响灰土强度的关键因素之一。胶结力学模型研究主要分为模型试验与数值模拟2个方向。模型试验研究即通过研发模型装置模拟土颗粒间的胶结作用进行试验研究。Delenne等[14]将由环氧树脂胶结的2根铝棒进行拉伸、压缩和剪切试验，得到简单应力路径下的强度参数，实现了胶结颗粒细观作用的宏观试验研究。蒋明镜等[15]在Delenne等[14]试验的基础上考虑法向作用力、胶结厚度及宽度等因素的影响，建立了三维微观胶结模型。Brockel等[16]和Kirsch等[17]研发了一套测试毫米级颗粒粒间固体胶结物荷载试验的装置，用特定的人工环境来对自然环境下晶体或盐类颗粒的形成与生长进行描述。数值模拟研究方面，Shen等[18]采用离散元法对2个大颗粒和若干个小颗粒组成的胶结球体进行了不同加载方式下的数值仿真研究。蒋明镜等[19]采用离散单元法对胶结砂土模拟直接拉伸、单双轴压缩及巴西试验。刘芳等[20]采用有限元方法模拟了圆形胶结颗粒微观抗压试验，分析了粒间胶结厚度和宽度对粒间法向接触应力分布和抗压强度的影响。

上述成果对研究土体宏观力学响应与微观力学行为的关联性具有重要价值。鉴于灰土是由多相复合材料组成，石灰在水分的作用下与各相材料发生强烈的相互作用，对土体起到胶结固化作用，导致灰土的颗粒分布具有复杂性和随机性[21-22]，含水率的不同也会影响灰土胶结特性。简单的力学模型并不适用于灰土古建筑。以灰土为主要材料的土遗址、城墙等古建筑属于不可再生文物，研究灰土力学性能衰减机理不能以牺牲文物为代价。因此，以电镜图片为依据，建立符合灰土古建筑微观结构的力学模型具有重要意义。

本文以国家级文物开封城墙为研究对象，根据300倍SEM图像，将城墙灰土分为颗粒、胶结物和孔隙三部分，基于灰土微观结构和图像处理技术，构建灰土分析模型。以颗粒粒径为依据，将粗颗粒定义为骨架颗粒，细颗粒定义为胶结物，通过试验方法建立胶结物的塑性损伤本构模型，以验证本文灰土模型的有效性，以期为研究以灰土为主要材料的土遗址、城墙等古建筑的力学特性提供参考依据。

1灰土力学模型的构建

土遗址、城墙等古建筑属于不可再生文物，取样受到一定限制，模型不能做得太大。通过建立相似模型模拟遗址实际工况，可以通过将实际问题缩小为一个模型来进行研究，同时将模型试验的结果反演到实际问题中去。相似模型试验的相似原理是指通过模型重现与原型相似的物理现象，即要求模型材料、性质及几何形状等均遵循一定的规律[23-24]。土体的宏观力学性质是由其自身的微观结构所决定的，取少量遗址土即可获得多张电镜图片，依据电镜图片进行微观建模，对土遗址进行力学特性分析。

1.1电镜图像倍数选取

选择合适的电镜图像倍数对建模的准确性影响较大。图1为开封城墙土不同放大倍数的电镜图像，分别为300倍、500倍、1 000倍、2 000倍。由图1可知，土样在放大300倍的视野范围内可以清晰看到颗粒、胶结物和孔隙三相结构、颗粒形状、尺寸和分布状况；将倍数提高至500倍～2 000倍时，颗粒数目较少，不具有代表性，失去对微观尺度下的指导意义。由于每张电镜图像的颗粒、胶结物和孔隙分布不均，为避免拍摄电镜图像的随机性影响试验结果，在进行300倍下电镜图像扫描时，需根据“从左到右”及“从上到下”的原则对土样进行拍摄。

图1开封城墙土SEM图像Fig.1SEM Image of Kaifeng City Wall Soil

1.2电镜图像处理

1.2.1 颗粒划分及提取

通过对土样的电镜图片进行处理，对颗粒进行划分和提取。土颗粒是具有连续且相近的阈值，成颗粒状分布[25]。为使处理后的电镜图像与真实土样的结构保持一致，本文基于专业图像处理软件Image-pro plus 6.0 (IPP6.0)、Photoshop、Vmde软件对300倍的SEM电镜图片进行处理和分析。颗粒图像处理与分割流程如图2所示。

图2颗粒图像处理与分割流程Fig.2Particle Image Processing and Segmentation Process

(1)通过IPP软件中的阈值工具设置阈值参数(90～105)作为颗粒区间，对图像依次进行空间刻度校准、图像预处理、二值化处理、数据输出[26]。

(2)使用Photoshop CC2018软件选择去杂化、锐化及加强图像工具，消除颗粒边缘与内部杂点，并做进一步细化处理，使用笔刷工具对照底图二维图片分别进行添加和删除手动操作，处理颗粒与胶结物黏连等情况。

(3)经Photoshop处理后的图片导入Vmde软件进行图像分割、编辑分割、平滑处理等操作，实现颗粒数据的分割提取，并导出dxf文件。ABAQUS软件可导入dxf文件进行建模。

1.2.2 土样孔隙率确定

通过对电镜图像进行分析，以灰度幅值的范围为依据，对组成像素进行分析，确定土样孔隙率。通过MATLAB软件数值分析功能，以SEM图像的阈值为研究对象，以积分为手段，通过三维层面对土体内部实际孔隙率进行求解[27]。颗粒及孔隙分布计算原理如图3所示，其中Y为白色灰度，Yi为黑色灰度。

图3颗粒及孔隙分布计算原理Fig.3Calculation Principle of Particle and Pore Distribution

通过灰度图像阈值为媒介，将SEM图像进行区域划分，积零为整求得所需数据。三维孔隙率V3D如式(1)所示。

(1)

式中：S为每一个像素的面积，选择单位像素，则S=1；M为图像灰度最大值，M=255；Di为第i个像素对应的灰度；N为图像像素的总个数。

运用MATLAB软件编写提取三维孔隙特征参数程序，将图3(a)的电镜图像进行处理，得到三维图像，如图4所示。

图4SEM图像三维二值化处理Fig.4Three Dimensional Binarization of SEM Image

1.3灰土模型建立

灰土模型是由颗粒、胶结物和孔隙组成的三相结构，根据以下规则对各相材料进行划分：

(1)采用IPP软件对SEM图像进行处理，通过阈值处理，将具有连续颗粒状灰色部分定义为颗粒，利用Photoshop、Vmde软件并提取颗粒边缘。

(2)由于颗粒间孔隙分布范围广，且形状极其不规则，定义存在困难。因此，在满足模型计算精度及效率的基础上对孔隙进行简化处理。颗粒间孔隙位置定义为3个颗粒形心连线的垂线交点处，形状为圆形[28]，孔隙大小范围根据文献[29]进行设置。若交点处为颗粒边缘和颗粒内部，则不进行孔隙定义，孔隙处理方法如图5所示，微观灰土的建模流程如图6所示。

图5孔隙位置处理方法Fig.5Treatment Method of Pore Position

图6建模流程Fig.6Modeling Flow

2模型本构关系及材料参数确定

将所建立的灰土模型导入ABAQUS有限元软件中计算，分别赋予模型中颗粒、胶结物、孔隙各自相应的材料属性。其中起支撑作用的颗粒以石英、长石颗粒为主，该类型颗粒硬度较大，施加一定的外荷载基本不会发生变形和破坏。有限元模拟中可以将刚硬的物体指定为刚体，不需要求解刚体的变形，可以减少计算时间，并提高计算效率。颗粒采用线性弹性体，材料参数取石英砂相应参数[30]。表1为颗粒与孔隙主要力学参数。胶结物采用塑性损伤本构模型(Concrete Damage Plasticity，CDP)模拟，该模型考虑了材料的抗压性能不同，因此需要确定胶结物相关力学性能参数，本文根据宏观胶结试验确定胶结物的应力-应变关系。

表1颗粒与孔隙主要力学参数Table 1Main Mechanical Parameters of Particles and Pores

2.1胶结试验

胶结物与土类似，因其结构的复杂性和不均匀性，需确定胶结物的力学性能参数。灰土可以看作由无数颗粒构成，在颗粒堆积的过程中，因为颗粒粒径不同，颗粒含量不同，颗粒构成矿物物质的不同，不可避免地对其形成的土体工程性质造成影响。较粗的颗粒如表面黏结石灰的砂粒主要承担骨架作用[31]；细小的粉粒、黏粒与部分石灰相互作用主要起胶结骨架颗粒及填充孔隙作用[32]。本工作以开封城墙灰土为研究对象，将粗颗粒定义为骨架颗粒，较细颗粒定义为胶结物，对灰土进行筛分试验，筛分粒径d分别为0.075、0.1 mm，制作含水率W分别为14%、15%、16%、17%、18%、19%的胶结材料土样进行压缩试验，建立含水率与胶结物强度特性的相关关系。

2.1.1 试验土样

试验所用土为开封城墙剥落土。通过对土样进行XRD矿物分析，确定该土样由石灰、沙子和黏土构成，且比例为3∶4∶3，采用了三七灰土，其基本物理参数如表2所示。

表2开封城墙土基本物理参数Table 2Basic Physical Parameters of Kaifeng City Wall Soil

2.1.2 胶结物粒径选择

开封城墙土颗粒级配状况是影响胶结物颗粒粒径选择的一个重要指标，也是划分土颗粒与胶结物类型的主要依据。开封城墙土的颗粒级配曲线如图7所示。

图7颗粒大小分布曲线Fig.7Particle Size Distribution Curve

从图8可以看出：粒径为0.005～0.075 mm的土样占灰土总质量的21.2%；粒径为0.005～0.1 mm的土样占比为28.9%。通过IPP软件处理得到的电镜图片进行颗粒粒径标注，如图8所示。由图8可知，颗粒粒径以大于0.075、0.1 mm为主。因此，本文将开封城墙土进行筛分，采用粒径小于0.075、0.1 mm的黏粒、粉粒与部分石灰作为胶结物进行分析。

图8颗粒粒径选取Fig.8Particle Size Selection

2.1.3 胶结物土样力学性能试验

将取回的开封城墙土使用密封袋保存，测定其天然含水率为14%。根据击实试验，确定土样最优含水率为19%。为保证试验土样一致性，将原状土过2 mm筛网备用。制备流程如图9所示。

图9土样制备流程Fig.9Soil Sample Preparation Process

(1)开封城墙土样制备：根据开封城墙土击实试验，确定土样的干密度。制备含水率为14%、15%、16%、17%、18%、19%的土样，干密度分别为1.47、1.49、1.50、1.52、1.55、1.57 g·cm-3，土样尺寸为φ39.1×80，体积为96 cm3，制备6组(每组4个平行土样)。

(2)胶结物含水率确定：考虑土体中含水率为均质的，胶结物的含水率选取为14%、15%、16%、17%、18%、19%。胶结物干密度确定：根据《土工试验方法标准》(GB/T 50123—2019)进行击实试验，在同一击实功下，确定土样在不同含水率下的干密度，制备土样并进行力学性能试验。

土样制备：通过筛分试验，将城墙土过0.075、0.1 mm筛网将胶结物分离出来，制备含水率W分别为14%、15%、16%、17%、18%、19%，0.075 mm土样的干密度分别为1.40、1.42、1.45、1.48、1.51、1.52 g·cm-3；0.1 mm土样的干密度分别为1.41、1.42、1.44、1.48、1.51、1.53 g·cm-3，土样尺寸为φ39.1×80，体积为96 cm3，分别制备6组(每组4个平行土样)。

利用单轴压缩试验测量胶结物与开封城墙土的抗压强度。采用无侧限压缩仪按照1 mm·min-1的轴向应变速率对开封城墙土及胶结物土样进行单轴压缩试验。

图10为开封城墙土与胶结物的应力-应变(σ-ε)关系(A为粒径d<0.075 mm，B为粒径d<0.1 mm)。由图10可知：在一定含水率范围内，开封城墙土的抗压强度大于胶结物抗压强度，两者的力学特性均为非线性；0.075、0.1 mm粒径胶结物的应力-应变变化规律基本保持一致。

图10土样与胶结物轴向应力-应变关系曲线Fig.10Relationship Between Axial Stress and Strain of Soil Sample and Cementitions

2.2胶结物塑性损伤模型

选用的ABAQUS塑性损伤模型(CDP模型)属于连续介质塑性损伤模型，是依据Lubliner等[31]、Abu-lebdeh等[32]提出的损伤塑性模型确定的，适用于模拟准脆性材料因为损伤而导致的不可逆的材料退化，考虑单向加载、循环加载及动态加载等情况，分析材料承受拉压能力的不同[33]。胶结物属于脆性材料[34-35]，其结构内部的特点是多样性的，存在刚度退化问题，也兼具塑性变形性质，产生的原因是各种细微孔洞、裂隙的延伸，在外应力作用下将产生细微损伤并逐渐扩大导致材料破坏[36]。

2.2.1 CDP模型塑性参数

CDP模型塑性参数由CDP统一硬化模型本构表示，统一硬化模型参数如表3所示。

2.2.2 CDP模型压缩损伤参数

CDP模型压缩损伤参数由单轴压缩下应力-应变关系确定。通常拟合岩土材料的应力-应变曲线，即寻找能在一定误差范围内描述该曲线的表达式，胶结物的力学性质表现为非线性，应力-应变关系复 杂。本文基于易达等[40]提出的岩土材料函数表达式，建立不同含水率的胶结物表达式，见式(2)；通过对d<0.075 mm、d<0.1 mm粒径的胶结物应力-应变曲线(图11)进行函数式拟合得到应力-应变拟合参数值，如表4所示。

表3统一硬化模型参数Table 3Parameters of Unified Hardening Model

图11胶结物应力-应变拟合曲线Fig.11Stress-strain Fitting Curves of Cementitions

(2)

表4不同含水率胶结物应力-应变拟合参数值Table 4Stress-strain Fitting Parameters of Cement with Different Water Contents

式中：εc为单轴压缩下峰值应力对应的应变；αa、αb、αc、c1和αd、αe、c2分别为应力-应变曲线上升段、下降段的相应拟合参数。

CDP模型中用来表示压缩损伤指标D的公式为

(3)

式中：E0为初始弹性模量。

3单轴压缩有限元数值模拟分析

3.1单轴压缩仿真模型

基于SEM图像建立模型，进行二维单轴压缩试验数值仿真，以验证所提出的灰土分析模型的有效性。通过对土样进行电镜试验，根据电镜图像在ABAQUS软件中进行建模。图12为含水率15%(S1～S4)、19%(S4～S8)土样的单轴压缩模型。

图12单轴压缩仿真模型Fig.12Simulation Model of Uniaxial Compression

3.2设置力学参数和网格划分

有限元模型采用ABAQUS/Explicit求解模块，划分的网络单元类型均采用CPE4R单元，该单元类型为四节点双线性平面应变四边形单元，减缩积分，使用该单元类型的好处是在涉及网格扭曲问题尤其是较大的网格变形时，可以使用细网格部分的线性减缩积分单元。有限元模型建立后，对模型的下边界进行位移约束。在进行数值模拟时，弹塑性损伤模型中的材料会出现非线性的性质，当荷载逐步接近极限承载力时，若给定荷载步过大，将会导致迭代发散，收敛困难。因此，将增量步大小的初始值设为0.001，最小值为1×10-12。

根据电镜图像尺寸大小，建立0.317 mm×0.647 mm的二维模型，将几何实体设置为可变性的实体。土体上方设置长0.317 mm、厚0.02 mm的刚性板。在定义相互作用时，板定义为接触的主面，土体上部设定为从面，法向行为中设置为硬接触，摩擦公式选为罚函数，摩擦因数取0.492。对土体单轴受压试验采用位移法进行加载，通过改变位移加载幅值Δ的大小来进行模拟，分别为0.003、0.006、0.009、0.012 mm。采用ABAQUS软件进行数值模拟时，对模型下边界施加X向、Y向及转角的约束。

3.3单轴压缩仿真模型验证

本文将单轴压缩仿真模型与室内单轴压缩试验进行对比分析。采用所建立的仿真模型模拟含水率为15%、19%的土样，将得到的S1～S4、S5～S8应力-应变曲线分别与宏观室内试验结果进行对比，结果如图13所示。可以发现，模拟与试验应力-应变曲线结果吻合良好，初步验证了所建立模型的可行性及参数选取的合理性。

图13模拟应力-应变曲线与试验结果对比Fig.13Comparison of Simulated Stress-strain Curves and Test Results

图14为微观数值仿真模型模拟破坏结果与宏观单轴试验结果的对比。从图14可以看出，模拟与试验破坏模式相同。在受力初期，试件裂缝开始萌芽，此时裂缝数量较少，破坏呈现出集中于主干式的单一主裂纹破坏；随着压力持续作用，试件裂缝不断加深及扩展，并贯穿试件，即沿试件斜向贯通，形成斜向连贯的断裂带，直至试件最终破坏。

图14模拟破坏结果与试验结果对比Fig.14Comparison of Failure Results Between Simulation and Test

3.4仿真模型压缩破坏

图15为试件S1在不同加载位移作用下的竖向位移云图。由图15可知，土体位移呈规律性变化，胶结物在压力作用下产生挤压变形。由于土体内部各相材料分布不均，土体位移呈不均匀变化，土体两侧出现“凹凸不平”的现象。在变形初期，位移主要集中在试件上方两侧；随着荷载的不断增加，土体的不均匀现象加剧，位移继续增大，局部位移发生突变，这说明内部各相材料的不均匀性导致土体内部应力分布不均，最终导致土体破坏程度不同，即土体的宏观力学性质是内部微观结构变化的外在宏观表现。

图15不同加载位移的竖向位移云图Fig.15Vertical Displacement Nephograms of Different Loading Displacements

图16不同加载位移的竖向应力云图Fig.16Vertical Stress Nephograms of Different Loading Displacements

图16为试件S1在不同加载位移作用下的竖向应力云图，其中，“+”为拉应力，“-”为压应力。分析表明，试件的应力整体分布不均，云图分布呈现竖向递进关系。竖向应力主要分布在试件两侧及孔隙处,中间部分应力较小，表明在单轴压缩荷载作用下土体试件内存在应力集中现象。

图17不同试件在不同位移作用下的损伤云图Fig.17Damage Nephograms of Different Specimens Under Different Displacements

图17为试件S1～S8的损伤示意图。由图17可知，土体的内部组分损伤(用损伤因子ω表征，ω=0表示无损伤，ω=1表示完全破坏)是相对不均匀的，土颗粒与胶结物界面区最先出现损伤，接着损伤以此向胶结物、孔隙继续发展，并在孔隙周围出现应力集中现象，影响损伤带的发展路径。随着荷载的持续增加，损伤区域不断扩展和演化，损伤集中在土体两侧，最终形成斜向连贯的塑性损伤区域，塑性损伤带交叉存在，土体整体产生软化行为直至破坏。

4结语

(1)以颗粒粒径为依据，将粗颗粒定义为骨架颗粒，细颗粒定义为胶结物，通过室内胶结试验，在一定含水率范围内， 0.075、0.1 mm粒径胶结物的应力-应变变化规律基本保持一致，建立了胶结物与含水率的函数关系式。

(2)根据灰土的胶结特性，将土体划分为颗粒、胶结物及孔隙组成的三相结构。通过300倍SEM图像重建法确定了颗粒的形状、大小及分布状况，建立了可以反映真实土体的微观有限元模型。

(3)通过对灰土力学模型进行单轴力学性能仿真，发现试件在单轴压缩荷载作用下，颗粒与胶结物界面处及孔隙周围首先出现损伤破坏并产生裂缝，其主要原因是该区域存在应力集中现象；随着荷载继续增大，损伤区域向外扩展延申，形成斜向连贯的断裂带，最终导致试件破坏。

(4)通过微观结构仿真模型模拟单轴压缩试验与宏观室内试验进行对比，验证了所建立力学模型的真实有效性，为研究以灰土为主要材料的土遗址、城墙等古建筑的力学特性提供参考依据。