微观尺度下州桥遗址土干缩劣化机理分析

岳建伟，黄轩嘉，赵丽敏*，王自法

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

受赋存环境的影响，自然界中土体的干湿状态总是在不断变化，土受组成材料（如颗粒组成、矿物组成和骨架结构）的限制及外部环境的影响，易造成失水收缩，力学性能降低。黄河中下游土遗址组成以粉质黏土为主，受外界环境反复作用，其表层裂缝不断扩展并逐渐贯通，这些裂隙不仅会对土遗址的物理特性产生影响，还引起其内部微观颗粒约束机制的失效。

土作为准脆性多相材料，由于其宏观试验步骤复杂、周期长、工作环境差和结构复杂，开裂过程的细节不能总是被观察和全面评估。为更好地了解土的失效机理，并从随机相位分布中评估其不确定性和可靠性，越来越多的学者重视微观尺度下土的损伤与断裂机理。

随着国内外学者逐渐将干缩裂隙发育的机理作为研究重点，反映其主要劣化机理的数值模型也得到长足发展，其中，有限元法对散体材料及破坏过程的模拟适应性较好，在模拟非饱和土方面已取得显著成果。程昊、崔贤等应用隐式积分算法开发了非饱和土扩展剑桥弹塑性模型程序，并在ABAQUS中对土体劣化及其稳定性进行模拟，取得了良好的效果。崔凯、陈毅、杨龙等研究发现土质材料粒度组分及微观结构是病害发育的基础，环境作用是病害发育的主要原因。Oualmakran、黄启迪、田俊峰等研究发现持续的失水过程导致土样内固-液交界处的弯液面曲率半径减小，使得作用在土颗粒骨架上的吸应力增大，在宏观上表现为拉力作用，导致了土体的收缩和开裂。韦秉旭等运用计算机图像处理，对不同干湿循环次数影响下土样表面裂隙特征进行了描述和定量分析，探讨干湿循环次数对土样开裂规律的影响。朱磊等通过土壤裂缝演化试验提取裂缝形态特征，并基于胡克定律建立了具有较高精度的2维土壤开裂模型。

上述学者主要以试验及本构关系进行了土样的劣化分析，并未将试验结果与现象进行数值化归纳。针对土遗址等存在特殊性的研究对象，存在无法大量取样和无法进行大量的筛分试验等难题。本文利用课题组前期基于SEM图像建立的2维数值模型，建立多个干缩试验模型并进行裂隙演化模拟，使用玫瑰图和位移场矢量图分析多个方位的劣化情况，更精准模拟土样开裂位置情况，并在宏-微观尺度下分析归纳州桥遗址土劣化机理、演化规律及其影响因素。研究方法和结果为实时跟踪分析土遗址的干缩劣化进展提供了理论方法，也为进一步探索土遗址保护技术提供了依据。

1 干缩试验及微观力学参数

1.1 干缩试验

试验用土取自河南省开封市州桥遗址。按照土工试验要求，将取回的土样装入塑料袋密封保存。按照《土工试验方法标准》（GB/T50123—2019）进行土工试验，土样物理力学性质如表1所示。

表1 土样基本物理指标
Tab. 1 Basic physical parameters of soil

液限/% 塑限/% 塑性指数38.5 23.2 15.2微毛细孔隙率/% 天然干密度/(g·cm-3) 天然含水率/%41 1.67 12.8

为掌握州桥遗址土的裂隙发育特征，主要研究土样干缩裂隙发展规律。依据《土工试验方法标准》（GB/T50123—2019）中第26节收缩试验要求，试验流程如下：1）将土样放入烘箱烘8 h；2）在直径61.8 mm的模具中制作土样，并称取环刀加土样的质量；3）将土样放入30℃烘箱中进行脱水，直到土样质量不再改变，停止试验，此为第1次循环；4）使用喷壶对第1次循环后的试样喷水，加至环刀加土样的质量与初始重量一致停止；5）重复步骤3）、4），直至完成3次循环。试样及模具如图1所示，试验开裂结果图通过Image Pro Plus 6.0（IPP）软件进行二值化处理。

1.2 土样力学特性

州桥遗址土的颗粒级配曲线见图2。从图2可以看出，粒径为0.005～0.074 mm颗粒的质量占比为53%，粒径为0.074～2.000 mm颗粒的质量占比为40%，并且在使用实验室最小粒径的0.02 mm筛进行筛分时已经很难再筛分出土样。通过SmileView软件对州桥遗址土的电镜图片进行土颗粒粒径标注（图3）。综上，考虑到实际试验效率，本文假定粒径小于0.074 mm的颗粒为细颗粒，粒径大于0.074 mm为粗颗粒。

图3 州桥遗址土电镜图（SmileView软件处理后）Fig. 3 Electron microscope of soil at Zhouqiao site (After processing by smileview software)

通过对遗址土样进行筛分试验，筛分出粒径小于0.074 mm的细颗粒和粒径大于0.074 mm的粗颗粒，然后依据《土工试验方法标准》（GB/T50123—2019）第20节无侧限抗压强度试验要求进行试验，最终得到粗颗粒和细颗粒土样的应力应变曲线，见图4。

图4 粗粒土与细粒土应力-应变Fig. 4 Stress-strain curves of coarse-grained soil and fine-grained soil

1.3 颗粒间的收缩力

宏观尺度下，土样应力-应变关系是土体物理特性最大影响因素；微观尺度下，表面张力对于土样干缩裂隙的产生影响最大。李靖等通过研究发现，大于10 nm的土样颗粒间的水分主要以液桥的形式存在，且孔内的水分会对土样稳定性产生影响。因此，在保持土样宏观物理特性不变的情况下，将不同表面张力作为梯度影响因素，得出微毛细孔隙间力与含水率之间的关系，为数值模型提供理论依据。

图5为土颗粒受力简图。图5中：

r

和

r

为互反弯月面的两个曲率半径；θ为水与固体表面的接触角；

R

为土颗粒半径；

φ

为弯液面水膜在球形土颗粒上的最大覆盖范围角，即饱和角；α为表面张力系数；

F

为毛细作用力；

T

为湿吸力。

图5 土颗粒受力简图Fig. 5 Soil particle force diagram

相同大小的两个土颗粒之间填充了凹面水膜，通过观察可以发现，两个土颗粒之间的液膜凹凸面上的水环是不规则体， 该凹凸面所带来的毛细作用力会造成土颗粒之间产生孔隙间力

F

，孔隙间力

F

由基质吸力

F

和湿吸力

T

组成，即：

本文土颗粒的平均粒径为0.074 mm，在20 ℃下，纯水的表面张力系数α为7.284×10N/mm，接触角为10°。最终通过对式（1）进行拟合可得微毛细孔隙间力与含水率的变化，如图6所示。

图6 微毛细孔隙间力与含水率关系[23]Fig. 6 Relationship between micro capillary pore force and water content[23]

2 有限元模型的构建

本文利用课题组前期基于SEM图像建立的2维数值模型，以宏观试验结合数值模拟的方式对州桥遗址土干缩裂隙机理进行研究。

模型尺寸为500倍镜头下拍摄的SEM照片尺寸大小，即644 µm×431 µm。在建模的过程中，将遗址土微观结构划分为3相：粗颗粒、细颗粒和微毛细孔隙。

依据压汞试验得到的土样微毛细孔隙率找出最符合真实土样结构的阈值，采用IPP软件进行图像处理。将处理后的图片导入到AutoCAD软件进行各组分边界划分，最终在ABAQUS软件内进行建模，处理流程如图7所示。

图7 建模流程Fig. 7 Modeling process

土颗粒及颗粒间胶结材料依据课题组前期理论推导及试验结果进行属性定义。数值模型依据表面张力作用机理施加作用力（图8），选用面积力进行施加荷载，表面张力大小依据图7确定。模型网格单元为CPE4R，对模型下边界进行

X

、

Y

及转角位移约束。最终，本文通过图9中的4张电镜图片建立干缩试验数值模型，如图10所示。

图8 表面张力作用机理Fig. 8 Mechanism of surface tension

图9 电镜图Fig. 9 SEM images

图10 微观仿真数值模型Fig. 10 Micro simulation numerical models

3 试验结果及分析

3.1 干缩试验结果分析

图11～13分别为干湿循环1、2、3次的试验土样表面开裂结果。

图11 干湿循环1次土样表面开裂结果Fig. 11 Surface cracking results of soil samples after one dry wet cycle

由图11可知：试样初始状态的表面平整；放入烘箱1 h后，试样表面开始变粗糙，随后试样表面不断出现微细裂缝并逐渐延伸和变宽；紧接着，垂直初始裂缝的微细裂缝也开始相继出现，如图11（b）所示；最终裂隙网络发育成形，如图11（c）所示。通过对1次、2次、3次干湿循环裂缝进行对比，发现干湿循环次数增加，主裂缝之间的距离有所增加，裂缝的数量明显增加，裂缝的分布越来越均匀（图12、13）。土样开裂的原因是由于土样表面存在着大量的微孔洞和微裂隙，在反复干湿作用下，这些微缺陷开始萌生、扩展汇合，直至形成宏观裂隙；干湿过程中温度应力的反复作用，也会导致土样颗粒的表面张力发生变化，造成不均匀收缩与膨胀，引起土样内部的微裂隙扩展，最终表现为土样的结构损伤。

图12 干湿循环2次土样表面开裂结果Fig. 12 Surface cracking results of soil samples after second dry wet cycle

图13 干湿循环3次土样表面开裂结果Fig. 13 Surface cracking results of soil samples after third dry wet cycle

3.2 数值模拟干缩试验机理分析

3.2.1 损伤云图结果分析

图14为模型1不同阶段下损伤云图，图15为分析步时间为1.00 s时模型1～4的损伤云图与试验结果对比。图14、15中，通过将微观模拟结果拼接成宽1.3 mm、长2.0 mm的矩阵图，使模拟结果从微米级及提升到可筛分的毫米级，旨在与宏观试验结果进行对比分析。

图14 模型1不同阶段下损伤云图Fig. 14 Damage nephogram of model 1 at different stages

图15 模型1～4损伤云图与试验结果对比Fig. 15 Comparison of damage nephogram of model 1～4 with test results

图14中蓝色部分为细颗粒，不规则白色部分为粗颗粒，圆形白色部分为微毛细孔隙，红色部分为损伤裂隙。限于篇幅及工作的相似性，本文仅将模型结果（图14模型1不同阶段下损伤云图）与干缩试验结果（图11～13）进行对比分析。州桥遗址土裂隙演化过程可归纳为：

1）张拉应力场的形成。在干燥过程中，土颗粒间不规则微毛细孔隙分布的微毛细孔隙水会通过微毛细孔隙排出，使粒间微毛细孔隙产生表面张力，进而使微毛细孔隙表面产生收缩变形。当干燥过程到达一定程度时，会在粒间孔隙表面形成张拉应力场，如图14（b）所示。

2）单根主裂隙的形成。通过数值模拟（图14（c））可以发现，随着土样局部发生应力集中现象，当张拉应力的大小超过土颗粒间的抗拉强度时，土样表面薄弱点处会萌发裂隙。在宏观试验结果中，土样均先产生单根主裂隙，并随着土样脆弱面不断延伸，如图11（a）、12（a）和13（a）所示。

3）二级、三级裂隙的形成。通过数值模拟结果（图14（c））发现在第1条裂隙产生后，裂隙延伸端的应力集中使得脆弱裂隙面上产生新的拉应力场，当其路径上的土颗粒位移受到限制或者发生应力集中时，便会衍生新的裂隙，每一次新裂隙的产生都会使得表面累积拉应力得到释放。结合宏观试验结果（图11（b）、12（b）和13（b））发现，在此阶段裂隙发展速度最快，并且在次生裂隙与邻近主裂隙相交时会停止延伸，但次生裂隙的宽度会继续增加。当次生裂隙的宽度值接近主裂隙时，在裂隙的薄弱点处会继续衍生出三级裂隙。

4）裂隙网络的形成。结合数值模拟结果（图14（d））和宏观试验结果（图11（c）、12（c）和13（c））发现，新衍生的裂隙会向已生成的裂隙靠拢。随着主裂隙和二级、三级裂隙路径不断贴近，最终州桥遗址土会被分割成若干个区域，形成裂隙网。

应用IPP软件对模拟结果和试验结果的裂隙率进行统计发现：数值模型裂隙率结果分别为15%、18%、14%及13%，均值为15%，与宏观试验15%的裂隙率结果吻合；模拟结果和试验结果的裂隙发展规律也较为吻合，通过前文分析可知，数值模拟法和宏观试验法所得到的裂隙发展过程均为：单根主裂隙、二级和三级裂隙，至裂隙网。两项结果对比的一致性表明基于微观尺度进行土遗址的干缩劣化机理的建模和分析方法是可行的。

3.2.2 位移云图结果分析

图16为模型1在不同阶段下轴向位移云图。分析图16发现，土样位移变化较为规律，整体位移呈现分层变化，其中底部位移最小，上部位移较大。在干缩过程中，试样内部位移分布不均衡，结合宏观试验结果（图11、12和13）发现，二级和三级裂隙发育最快；在数值模拟分析中发现，在分析步时间0.52 s时局部位移发生突变，位移处于在2.498×10～6.246×10mm的部分最多（在图16中主要为天蓝色区域），比分析步时间0.37 s时1.249×10～2.498×10mm的位移变化幅度接近3倍，表明二级和三级裂隙发育速度要比张拉应力场的形成速度快。

图16 模型1不同阶段下轴向位移云图Fig. 16 Axial displacement nephogram of model 1 at different stages

图17为模型1～4轴向位移云图。分析图17发现，造成上述现象的原因是从微毛细孔隙水中逃逸出来的水分子需要穿越上层的干燥土样才能到达蒸发面，水分子的迁移路径因此而延长；此外，微毛细孔隙中气泡的存在还会干扰水分原有的迁移路径，增加水分迁移到蒸发面所需的时间，造成上部土样先于下部土样开裂。

图17 模型1～4轴向位移云图Fig. 17 Model 1～4 axial displacement nephogram

图18为土样内部位移矢量和位移趋势线，图19为土样内部损伤云图，图20为土样表面位移矢量和位移趋势线，图21为土样表面损伤云图。图19、21分别为图13中土样内部损伤云图与土样表面损伤云图；图18、20给出了图19、21中裂纹和新发展的裂纹（红色部分）附近的位移矢量。为了便于分析和解释损伤裂纹的趋势，在图18、20的图像中，分别绘制了一般位移趋势箭头，并与图19和21的损伤图像相对应。从位移矢量图的趋势线分析可知，裂纹萌生阶段几乎完全由微拉伸裂纹组成，而裂纹发展阶段开始由微拉伸裂纹和微剪切裂纹组成。通过对已有裂纹和新产生裂纹周围位移场分析，确定了两种位移场类型，分别为拉伸位移场与拉伸-剪切位移场；而拉伸位移场既与拉伸裂纹的发展有关，也与混合拉伸-剪切裂纹发展有关。从图18、20中可以看出，表面张力的方向并不固定，是随着裂纹的发展而变化的，而且逐渐由拉伸作用力向混合拉伸-剪切作用力发展，如图19、21所示，随着张拉场的形成，主裂缝受到混合拉伸-剪切作用力，会垂直于受力面向外延伸开裂。图18与图19中所呈现的结果与Zhang等的研究结果相近。

图18 土样内部位移矢量和位移趋势线Fig. 18 Displacement vector and displacement trend lines in soil

图19 土样内部损伤云图Fig. 19 Internal damage nephogram of soil

图20 土样表面位移矢量和位移趋势线Fig. 20 Displacement vector and displacement trend lines of soil surface

图21 土样表面损伤云图Fig. 21 Nephogram of soil surface damage

土样内部与土样表面位移矢量和位移趋势线的发展也出现了不同，对比图18、20发现，在裂纹萌生阶段，土样内部位移矢量和位移趋势线指向作用力中心，而土样表面位移矢量和位移趋势线指向下层土样。这说明土样内部受到表面张力作用，微毛细孔隙收缩，而导致上层土样发生向下的沉降趋势，裂隙也沿着颗粒表面发展。在裂纹发展阶段，土样内部逐渐由拉伸位移场向混合拉伸-剪切位移场过渡，主裂缝垂直于受力面向外延伸开裂；表层土样由于裂缝的产生打破了区域应力场的平衡，使裂缝边缘的土颗粒失去相应的约束，在外力作用下，重新排列到断裂面的两侧，直到在新位置达到应力平衡，这也是导致土样表面开裂的根本原因。上述仿真结果与唐朝生等所提出的土样龟裂机理示意图趋势相同。

3.3 干缩对水平和竖向裂隙率的影响

为保证数值模拟结果的有效性，将仿真模型结果与宏观试验进行归一化处理，以时间与最大时间的比值（

x

=

t

/

t

）为横坐标，以裂隙率与最大裂隙率的比值（

y

=裂隙率/裂隙率）为纵坐标，绘制出时间同仿真模型结果与宏观试验无量纲裂隙率特征曲线，如图22所示，其变化趋势基本保持一致。

图22 数值模拟与宏观试验裂隙率归一化结果对比Fig. 22 Comparison of normalized fracture rate between numerical simulation and macroscopic test

4个模型的横向裂隙率、竖向裂隙率和总裂隙率的结果如图23所示。初始阶段裂隙率增长较慢，说明张拉应力场形成期较长；由于单根主裂隙的形成使得裂隙率有发生突变的趋势，随后二级和三级裂隙形成，使裂隙率实现跳跃式增长，土样发生大规模不稳定性劣化，最终裂隙网络逐渐稳定，裂隙率不再增长。

图23 裂隙率数值模拟结果Fig. 23 Numerical simulation results of fracture rate

结合图15、17、23可以发现：在裂隙发育初期，竖向裂隙率的发展速率要比横向裂隙率快；上层土样水分蒸发快，而下层土样由于受到水分子的迁移路径长的原因，导致上部土样先于下部土样开裂，并造成土样会沿着水分的迁移路径开裂，给竖向裂隙的发展提供了条件。随着水分迁移路径趋于稳定，竖向裂隙的尺寸发展趋于缓慢，水平裂隙成为土样裂隙发展的主要因素。

将数值模型的中心设置为基点，每间隔45°划分为一部分，可将图15均分成8个部分。采用IPP软件计算图15中4个数值模型中每一个部分所含裂隙（黄色部分）在每个区域的占比，从而绘制出8个方向的裂隙率玫瑰图，如图24所示；裂隙率结果如表2所示。

表2 数值模拟裂隙率占比
Tab. 2 Proportion of numerical simulation fracture rates

%裂隙率位置 模型1 模型2 模型3 模型4 均值模型上半部分 45 46 44 45 45模型下半部分 55 54 56 55 55 45°～135°区域与225°～315°区域 55 54 52 46 52 315°～45°区域与135°～225°区域 45 46 48 54 48

图24 裂隙率数值模拟玫瑰图Fig. 24 Rose diagram of fracture rates in numerical simulation

由图24和表2发现，4个模型的上半部分和下半部分的裂隙率分别占总裂隙率的45%和55%，上半部分土样的裂隙率普遍比下半部分低，并且裂隙多集中于45°至135°区间内和225°至315°区间内。结合图15发现,孔隙较多的方位裂隙率也更大，例如图24（a）中90°至135°区间、图24（b）中225°至270°区间、图24（c）中270°至315°区间及图24（d）中0°至45°区间，结果表明微毛细孔隙是造成土样开裂的主要原因之一。

4 结 论

基于SEM图片建立州桥遗址土在微观尺度下的数值模型并进行裂隙演化模拟；使用玫瑰图和位移场矢量图分析多个方位的裂化情况，更精准地模拟土样开裂位置情况并减少分析时土样颗粒不均匀对试验结果带来的影响；结合宏观试验结果，在宏-微观尺度下分析归纳州桥遗址土劣化机理、演化规律及其影响因素，主要得到以下结论：

1）州桥遗址土裂隙演化过程可归纳为拉应力场的形成、单根主裂隙的形成、二级水平、三级竖向裂隙的形成、裂隙网络。拉应力场与单根主裂隙的形成阶段持续时间最长；二级水平裂隙与三级竖向裂隙的形成阶段持续时间最短，且裂隙长度快速发展；裂隙网络属于发育减缓趋稳阶段。

2）裂隙类型产生的机制可归纳为局部剪切、局部拉伸及局部混合拉伸-剪切3种形式；随着裂纹的发展，逐渐由拉伸作用力向混合拉伸-剪切作用力发展，主裂缝垂直于受力面向外延伸开裂。

3）由于不同深处土样水分子的迁移路径长度不同，上层土样较更深处的土样更易产生较大的收缩变形，并给竖向裂隙的发展提供了条件。随着水分迁移路径趋于稳定，竖向裂隙的尺寸发展趋于缓慢，水平裂隙成为土样裂隙发展的主要因素。

4）分析裂隙率玫瑰图发现，上半部分土样的裂隙率普遍比下半部分低，并且裂隙多集中于90°与270°方向上，表明微毛细孔隙及土样内部收缩过大是产生裂隙的原因。

5）数值模拟与宏观试验裂隙率归一化结果对比的一致性表明，基于微观尺度进行土遗址的干缩劣化机理的建模和分析方法是可行的，本文的研究方法为进一步研究土遗址的保护技术奠定了基础。