基于有限元遗址土应力-应变关系的方法研究1)

岳建伟 黄轩嘉 赵丽敏 孔庆梅 陈 颖 王自法

*(河南大学土木建筑学院，河南开封 475004)

†(开封市不可移动文物修复与安全评价重点实验室，河南开封 475004)

土遗址作为古代建筑物的遗留部分，受组成材料的限制及外部因素的影响，易造成失水收缩，力学性能降低，最终产生劣化现象[1]。这些裂隙不仅会对土遗址的工程特性产生影响[2]，还会对遗址的稳定性及长期保存构成严重威胁，这也是土遗址保护领域公认的工程“疑难杂症”。

学者们对遗址土的宏观力学性能进行了一定的研究[3]，而细观尺度下的物理力学性能研究成果较少，遗址土的劣化机理缺乏系统的理论基础支撑。复杂力学行为数值仿真研究方向的飞速发展，能够有效避免物理试验中人为及环境因素的影响，满足工程领域研究遗址土基础特性的目的，而建立遗址土组成材料的本构关系及数值模型是满足工程应用的首要工作，也是影响数值模拟结果的重要因素。

岩土材料的力学性能极易受工程环境的影响，由此产生的裂缝显著影响土壤的力学特性[4]。冯永等[5]基于室内试验及有限元数值模拟研究了土的细观多相变形机理及其宏、细观变形特性之间的相互关系；郭梦圆等[6]通过对土进行单轴抗压强度实验，并对S–M 模型中各参数的意义进行修正，得到土在不同温度及载荷加载等级下的应力-应变曲线和改进的S–M 蠕变显式模型。徐安全等[7]通过压缩试验研究了压缩模量的变化规律并利用MATLAB 三维拟合得出压缩模量的表达式，更直观表现压缩模量的变化规律。

上述研究成果在一定程度上推动了宏、细观土样力学特性研究的发展，但鉴于土样自身的复杂性及多变性[8]，建立模型时，均对颗粒间的胶结物进行了一定的简化，同时也简化了不同环境作用下土样颗粒间胶结物的性能变化情况，这会使土样和土样各组成成分的应力变化规律存在偏差[9]。由于不能在既有土遗址本体上进行大量的取样以及进行反复的试验测试，本工作基于材料类比思想和广义相似理论，对比分析遗址土细观尺度和混凝土宏观尺度组成的特点，获得材料宏观和细观间的内在联系，推导颗粒间胶结物的计算公式，依据遗址土扫描电子显微镜(scanning electron microscope， SEM) 图像及少量的力学测试结果，以探索细观尺度下遗址土颗粒间胶结物的力学性能的确定方法。

1 研究思路

细观尺度下(10-6～10-9m) 的遗址土由土颗粒、颗粒间的胶结物、气泡三部分组成，而宏观尺度下(10-2m) 混凝土也由石子、石子间的砂浆、气泡三部分组成。石子和土颗粒为弹性体，受力过程中假设为处于弹性状态，相对于砂浆而言，土颗粒间胶结物的力学性能不清楚，本文以混凝土、水泥砂浆及遗址土的力学性能参数为依据，基于材料特性类比思想，根据相似理论及Hooke 定律推导出胶结材料力学性能参数[10-11]，为研究遗址土的力学性能和劣化奠定基础。

试验用土取自河南省开封市州桥遗址，按照《土工试验方法标准》(GB/T 50123–2019) 进行土工试验，土样基本物理指标液限为32.8%，塑限为18.7%，塑性指数为14.1，孔隙率为41%，天然干密度为1.67 g/cm3，天然含水率为12.3%，不同含水率下遗址土和混凝土的应力- 应变曲线如图1 和图2所示。

图1 遗址土应力-应变曲线

图2 混凝土应力-应变曲线[12]

混凝土与遗址土两者的力学性能满足相似理论的第一定理，主要有如下相似特征：

(1)分析图1 和图2 发现，含水率升高对两类材料的力学性能均有影响，并呈现出含水率越大，应力及应变值越小的情况。

(2) 图3 为宏观尺度混凝土材料示意图(10-2m)，图4 为细观尺度遗址土材料示意图(10-6～10-9m)，其中，ΔL为对应材料沿F方向的单位变形量。分析图 3 和图 4 发现，宏观尺度下混凝土整体性能是粗骨料及水泥砂浆共同作用的结果[13]，在细观尺度下遗址土整体性能是土颗粒及胶结物共同作用的结果[14]，两类材料在组分划分和构成规律、结构受力特点、力学特性等方面均较为相似。

图3 宏观尺度混凝土示意图

图4 细观尺度遗址土示意图

(3) 虽然两类材料的力学性质在数值上存在差距，但是两者的应力-应变历程规律较为相似，均属于准脆性材料，因为这类材料在载荷作用下，当自身应力、应变和弹性模量超过某一临界值后，材料会出现明显软化现象[15]，在实际工程中，对于遗址土和混凝土两类材料，当本体出现微裂缝后，其力学特性和工程质量会发生改变，为了使两类材料的力学性能发展规律更为接近，因此本文主要针对材料发生破坏前的性能进行分析，不再考虑后期的材料软化阶段，也满足实际工程需求。

(4)遗址土和混凝土结构在受力破坏后具有分形特性[16-17]。

2 模型参数的确定

基于遗址土细观尺度和混凝土宏观尺度组成的特点，通过材料类比思想和广义相似理论，推导颗粒间胶结物的计算公式，依据遗址土SEM 图像及少量的力学测试结果，通过下列方法探索细观尺度下遗址土颗粒间胶结物的力学性能的确定方法。

混凝土宏观尺度平衡方程与遗址土细观尺度平衡方程可表示为

式中，σx与σy分别为材料在X和Y方向的正应力分量，τxy为材料在X和Y方向的剪切应力分量，X与Y分别为材料在X和Y方向的体力。

两类材料的宏观力学与自身的材料组成密切相关，表1 为各物理力学参量相似常数，其中，Cl1为混凝土与遗址土几何长度l的相似常数，Cσ1为混凝土与遗址土应力σ的相似常数，Cε1为混凝土与遗址土应变ε的相似常数，Cγ1为混凝土与遗址土容重γ的相似常数。基于广义相似理论将表1 中相关式代入式(1)，可得

表1 各物理力学参量相似常数

并且为了让遗址土能反映出混凝土的应力状态，要求式(1)与式(2)平衡方程必须一致，并且强度极限与应力的量纲一致，因此混凝土与遗址土强度相似模数为

同理，可得水泥砂浆与胶结材料的平衡方程为

式中，Cσ2为水泥砂浆与胶结物应力相似比，Cγ2为水泥砂浆与胶结物容重相似比，Cl2为水泥砂浆与胶结物几何长度相似比，σ水泥砂浆为水泥砂浆材料应力，γ水泥砂浆为水泥砂浆材料容重，l水泥砂浆为水泥砂浆材料体长，σ胶结物为胶结物材料应力，γ胶结物为胶结物材料容重，l胶结物为胶结物材料体长，K2为水泥砂浆与胶结物相似模数。

基于相似理论第一定理可得K1=K2，则式(3)与式(4) 可联立为

对式(5) 进行简化变形可得遗址土颗粒间胶结物的应力关系式为

基于Hooke 定律可得材料中的应力与应变之间的关系为

式中，F为作用力，A为受力面积，ΔL为图3 所示各相材料沿F方向的单位变形量值，L为各相材料沿F方向的总变形量值。

由于粗骨料比水泥砂浆的弹性模量大了近一个量级，导致混凝土的宏观变形主要体现在粗骨料之间水泥砂浆接触状态的变化，可将混凝土材料沿F方向的单位变形量值ΔL混凝土视为由粗骨料及水泥砂浆材料两部分的沿F方向的单位变形量值ΔL构成，即

联立式(9)与式(10)并同时除以沿F方向的总变形量值L可得水泥砂浆与混凝土材料的应变比为

针对本工作的试验对象具有不可重复性，不能在既有本体上进行大量的取样，考虑借助于SEM 及CT 图像等细观结构图像结果进行研究(图5)。所以对式(11) 中沿F方向的总变形量值L取电镜图片单位宽度得到

图5 混凝土图像处理

式中，idem为混凝土中水泥砂浆占混凝土CT 图片的单位宽度面积的倒数。当分段数为1，2，···，10 时，水泥砂浆与混凝土电镜图片单位宽度面积比分别为0.35，0.37，0.48，0.51，0.49，0.42，0.57，0.41，0.47，0.53。可知idem= 0.46。

将宏观试验结果[19]中混凝土与水泥砂浆的应变比(图6)与利用式(12)推导所得到的应变比进行对比，发现理论与试验结果存在偏差。为修正混凝土与水泥砂浆应变关系存在的偏差，将宏观试验结果与理论推导结果的比值α定义为

图6 混凝土与水泥砂浆宏观试验应变关系

式中，εs，p为混凝土宏观试验应变值，εs，m为水泥砂浆宏观试验应变值，εl，p为混凝土理论应变值，εl，m为水泥砂浆理论应变值。

在土力学中，一般认为土颗粒是刚性的，土的宏观变形和破坏主要体现在土颗粒之间胶结物的接触状态的变化，则同理式(12)，胶结物与土样的应变关系为

式中，idem为遗址土中胶结物占遗址土电镜图片的单位宽度面积的倒数，由遗址土图像(图7)可知，当分段数为1， 2，···， 10 时，胶结构与遗址土电镜图片单位宽度面积比分别为0.64， 0.67， 0.57，0.74，0.78，0.65，0.61，0.51，0.63，0.53，可得idem=0.63。

图7 遗址土图像处理

依据式(13) 及式(14) 可得遗址土颗粒间胶结物的应变关系式为

将遗址土基本物理指标：遗址土样重度γ土=18.57 kN/m3，土颗粒重度γ颗粒=24.31 kN/m3，混凝土重度γ混凝土= 24 kN/m3，水泥砂浆重度γ水泥砂浆=20 kN/m3分别代入式(16)、式(6)及式(15)可求得遗址土颗粒间胶结物应力-应变曲线(图8)。

图8 遗址土与胶结物应力-应变曲线

3 细观尺度仿真模型的构建

土遗址大多处于地上属于无围压状态，并且土遗址的劣化大多集中于土体外表面，加之准脆性材料(concrete damage plasticity， CDP) 本构的计算结果不存在宏细观尺度上的差异[20]，所以本文基于SEM图像，并结合遗址土与胶结物应力-应变曲线进行二维单轴压缩试验数值模拟，以验证本文提出的遗址土颗粒间胶结物理论模型的有效性。

3.1 建模依据

图9 为宏细观尺度土颗粒结构分布，由图9 可知500 倍电镜图片土体颗粒清晰，不规则裂缝较为明显，在该尺度下可以观测到直径约5～25 μm 的土颗粒、孔隙的分布及无序排列的黏性颗粒；当提高电镜倍数到2000～5000 倍时，颗粒形状清晰，但颗粒数目较少，不具有普遍性；倍数提高到10 000～20 000倍时细观结构中颗粒形状模糊。综上，本工作依据500 倍状态下的电镜图建立遗址土的细观数值模型。

图9 宏细观尺度土颗粒结构分布

为避免拍摄电镜图片的随机性影响研究结果，在拍摄SEM 图片时，在满足筛分试验级别的基础上，对1 mm×1 mm 的面积域按从左至右及从下到上的原则进行取景，依据500X 下镜头拍摄的照片像素尺寸面积可得需拍摄电镜图片张数为4 张。

对图像进行二值化处理的关键是阈值的选择与确定，由于各种图像上目标物性质的不同以及灰度变化的多样性，使得传统的二值化方法难以取得理想的处理效果。为了使处理后的图像与真实土样结构吻合，应用Image-Pro Plus 软件对500X 放大倍数下的SEM 照片进行处理，为缩小单一阈值下计算结果随机性带来的误差，通过调整阈值大小所代表的切平面高低，以压汞试验孔隙率结果为依据确定黑白色度的面积比，找到最能真实反应土体细观结构的最佳阈值[21]，处理方式如图10，图中填充部分为土体颗粒，未填充部分为胶结物及孔隙，最终处理结果如图11 所示。

图10 不同阈值下的切片处理

图11 图像处理结果

将图像软件处理后的图片导入到AutoCAD 软件并依据以下规则对图11 进行各相材料边界划分：

(1)将灰色部分视为土颗粒，红色部分视为胶结物作为遗址土数值模拟建模宽泛规则，并忽略距离土颗粒形心较近的红色散点。

(2)在满足计算精度及效率的前提上，粒间孔隙中点位置为各个土体颗粒的形心连线的垂线交点[22]，形状为圆形[23]，孔隙大小范围依据文献[24] 进行设置。若存在交点处于颗粒边缘或颗粒内部，则不进行孔隙位置定义，孔隙位置具体处理方法如图12所示。

图12 孔隙位置处理方法

在建模时不考虑土体颗粒表面微孔隙的吸水作用，主要有以下两个原因：由于土体颗粒表面孔隙存在孔隙内气压高，颗粒间气压低的现象，并根据DLVO (Darjaguin Landau Verwey Overbeek) 理论，颗粒与水分子之间存在电荷层的相互作用，颗粒周围的水分子受细观力学作用仅在表面形成吸附水[25]。

3.2 模型建立

本文根据500X 电镜图尺寸大小，通过ABAQUS 软件构建长为0.644 mm 的刚性板和0.644 mm× 0.431 mm 的二维几何模型，以刚性板模拟试验中的压力设备，将几何体设置为实体，设置平面厚度为0.020 mm 并选取二维实体单元赋予截面属性，刚性板的弹性模量为206 GPa，泊松比为0.3。

数值模型由土颗粒、孔隙和胶结物三相组成，将土体颗粒定义为弹性刚体[26]，颗粒间胶结材料采用混凝土塑性损伤(concrete damaged plasticity，CDP)材料本构进行属性定义。该数值仿真模型中各相材料参数依据具体工况数值进行标定，土颗粒：弹性模量为2.0 GPa，泊松比为0.3；胶结物：弹性模量为21 MPa，泊松比为0.3。统一硬化模型参数：膨胀角Ψ = 30°，流动势偏移量ϵ= 0.1，双轴极限抗压强度与单轴受压极限强度之比fb0/fc0= 1.16，拉伸子午面上与压缩子午面上的第二应力不变量之比K=0.666 7，黏性系数μ黏=0.000 5。

本文采用位移加载法，对上压板顶部分别施加0.006 mm，0.009 mm 和0.012 mm 的竖直位移并采用平滑分析步进行加载。为确保土体与上压板能够真实地模拟接触面之间的滑移以及破坏等现象，在定义相互作用时将刚性板下边缘定义为接触主面，土样模型上边缘设定为接触从面，在法向行为中选为硬接触，在切向行为中摩擦公式选为Penalty (罚函数)，摩擦系数μ摩取0.492[27]。有限元模型建立之后，将模型下边界的X方向、Y方向及转角设置为固定端约束，上压板的X方向及转角设置为固定端约束。模型网格单元类型采用四节点双线性平面应变四边形单元，单元库选择动力学、单元族选择平面应变，网格全局尺寸为0.005。为防止给定载荷步过大，导致迭代发散，收敛困难，采用显示动力学进行计算。图13 为建模处理流程，图14(b)(d)(f)(h) 分别为含水率10%，12%，14%，16%土样的单轴压缩仿真模型结果。

图13 建模流程

图14 单轴压缩试验仿真模型

4 单轴压缩试验仿真验证及结果分析

由于篇幅有限以及工作的相似性，本文仅对模型1 进行不同阶段下轴向位移云图的分析。图15 为模型1 在不同竖直压缩位移下的轴向位移云图，图15 中图例“+” 为正方向，“-” 为反方向，图15 及图16 中参数量纲为mm。由图15 可得，土体位移变化比较有规律，整体位移呈现出分层依次变化，其中底部位移最小，而上部位移较大，因模型底部边界是固支约束，所以底部轴向载荷位移必为0。试块在单轴抗压过程中，其内部位移分布不均衡，且这种位移的不均匀性，随载荷的增加而加剧，在压缩位移持续作用下，局部位移会发生突变。图16 为模型1～模型4 在竖直压缩位移12 μm 下的轴向位移云图，通过对图16 中4 个模型轴向位移云图的分析发现，造成上述现象的原因是土体颗粒分布具有不均匀性，当有较大压缩位移时，土体颗粒与胶结物之间会发生搓动，土体内部产生较大位移，试样内部膨胀并造成内部应力过大，宏观表现为靠近中部的土体最先发生脱落，仿真模型结果与宏观试验规律表现一致。

图15 模型1 在不同竖直压缩位移下的轴向位移云图

图16 模型1 ～模型4 在竖直压缩位移12 μm 下的轴向位移云图

图17 为模型1 在不同竖直压缩位移下的轴向应力云图，图17 中“+” 为拉应力，“-” 为压应力，图17 及图18 中参数量纲为Pa。由图17 可知，试件压应力较大部分主要分布在中间部分，试件两侧压应力较小，整体应力呈现向外侧扩展的趋势，这也是试样中部发生较大位移的原因。此外，模型两侧的应力变化并不均匀，云图分布特征呈现出竖条递进关系，此外试块中部及两侧应力数值正负号不同，表示中间的应力最大，因为与边界区域相比，土壤是封闭的，在压缩过程中土体会发生横向膨胀的现象。图18 为模型1～模型4 在竖直压缩位移12 μm 下的轴向应力云图，通过归纳图18 中四个模型的轴向应力云图，发现在受压过程中，内部位移的不均匀性是导致试块内部应力复杂的主要原因，其在宏观表现为试件内部的破坏程度不同。

图17 模型1 在不同竖直压缩位移下的轴向应力云图

图18 模型1 ～模型4 在竖直压缩位移12 μm 下的轴向应力云图

图19 为模型1 在不同竖直压缩位移下的损伤云图，通过对图19 分析可得，试件损伤开始萌生并起始于试件底部两端，此时土体内部粘聚力让试块变形可恢复，试块表面无裂纹产生；随着压缩载荷的不断增大，土体内部基体产生少量微裂纹，产生不可恢复变形；随着损伤微裂纹进一步扩展、聚集，微裂纹不再是单一地某一条发展，而是沿着土体颗粒边缘全断面同时刻一起扩展。当应力达到一定值时，某些区域损伤程度开始加大，并且逐渐贯通，贯穿面外围出现主裂缝，最后形成几条贯通于整个试件的裂纹。图20 为模型1～模型4 在竖直压缩位移12 μm下的损伤云图，从图19 和图20 中可以看出，孔隙对土体的破坏有一定的影响，即孔隙周围容易出现应力集中现象导致损伤发生，而且孔洞对损伤路径的发展有影响，并影响试件最终的破坏面。

图19 模型1 在不同竖直压缩位移下的损伤云图

图20 模型1 ～模型4 在竖直压缩位移12 μm 下的损伤云图

图21 为宏观试验与仿真结果，图中图例括号内数字为土样含水率。本文通过对比理论模型与单轴压缩数值模型结果(图21) 发现数值模拟与试验应力-应变曲线结果吻合度高，初步验证了所建立模型的可行性及参数选取的合理性。所求得的结果已满足实际工程的需求，为下一步分析环境对土遗址劣化作用和土遗址安全评价奠定了基础。

图21 宏观试验与仿真结果

5 结论

本文基于材料类比思想、广义相似理论推导了颗粒间胶结物的计算公式，探索了细观尺度下遗址土颗粒间胶结物的力学性能的确定方法；基于遗址土500X 电镜图和图像重建法，构建了颗粒、胶结物和孔隙三部分组成的细观模型，采用位移加载方法获得了不同含水率遗址土的应力-应变关系，并与宏观试验结果进行对比，从细观角度对胶结颗粒力学特性进行了分析。主要结论如下：

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

(3) 通过对细观有限元模型进行单轴压缩试验，发现颗粒与胶结物接触位置与孔隙周围首先出现损伤，其主要原因是该区域存在应力集中现象；随着载荷增大，损伤区域向外扩展延伸，最终导致试件破坏。通过对比发现数值分析结果与宏观试验应力-应变曲线结果具有很好的一致性，表明本工作关于遗址土的研究思路和研究方法是可行的，也为取样困难的室内试验提供了新方法。