考虑土体结构劣化的高原山区桩-土体系数值分析

2022-08-19寇海磊侯王相荆皓陈琦李恒

湖南大学学报(自然科学版) 2022年7期

寇海磊，侯王相，荆皓，陈琦，李恒

（1.中国海洋大学工程学院，山东青岛 266100；2.国网山东省电力公司德州市陵城区供电公司，山东德州 253500；3.青岛市建筑工程管理服务中心，山东青岛 266071）

随着川藏铁路建设的持续推进，西部高原山区基础设施建设持续开展.高原山区桩基工程在服役过程中常年反复的冻融循环导致桩-土界面作用弱化，进而造成上部构筑物产生不均匀沉降、结构断裂等工程灾害，对上部结构物的稳定性带来了极大的挑战［1-4］.

冻融循环作用下寒区土体的物理力学性质变化始终备受国内外学者的关注［5-8］.Chamberlain 等［9］通过室内试验对冻融条件下寒区土体的渗透性以及微观结构进行了研究，发现土体在冻融过程中孔隙增大，且土体的原生结构遭到破坏；Hotineanu 等［10］通过室内模型试验探究了冻融条件下不同种类寒区土体力学性质变化规律，试验结果表明，冻融作用对土体的黏聚力影响较大，是影响土体稳定性的主要因素；苏永奇等［11］探究了冻融循环作用对青藏粉质黏土的力学参数影响，认为土体在冻融过程中发生了结构损伤，最终导致寒区土体力学性能产生不可逆的变化；赵茜等［12］通过室内三轴固结渗透试验探究了冻融循环条件下土体渗透性能的变化规律，结果表明土体的渗透系数随冻融循环次数的增加不断增大.虽然国内外学者对冻融条件下不同土体的物理力学性质变化进行了大量研究，但关于土体结构损伤对桩-土体系的影响及土体结构损伤模型的建立研究甚少.

国内外学者对高原山区桩-土体系相互作用进行了大量的物理模型试验［13-15］.因数值手段具有便捷性，数值模拟一直受各国学者的青睐［16-19］.Lai等［20］通过数值耦合模型对服役期内桩基冻胀过程进行了非线性分析；Lu 等［21］提出了一种预测冻土中单桩冻胀效应的数学模型，阐明了土体冻结深度对桩身轴力以及竖向位移的影响，但其并未考虑土体结构劣化对桩土界面相互作用的影响；Aldaeef和Rayhani［22-23］研究了桩基础在冻结、非冻结黏土中的荷载传递规律，并提出了一种基于冻结温度以及黏聚力因子预测桩基承载力的方法，但该方法未考虑桩-土界面相互作用随冻融循环不断弱化的因素；邢爽等［24］基于非线性计算模型对地震作用下季冻区冻土-结构动力相互作用进行了研究，但未考虑土体结构随冻融作用的劣化特性；唐丽云等［25］探究了不同温度条件下地下水位对多年冻土区桩基承载性能的影响，但未考虑常年冻融作用对桩基承载特性的影响；陈坤等［17］对多年冻土地区灌注桩桩身温度分布特性以及桩-土体系导热过程进行了研究，但未考虑土体因温度变化发生的冻胀融沉特性.综上，已有研究虽对考虑多场耦合条件下桩-土体系承载特性变化规律进行了研究，但对冻融作用下土体结构劣化所导致的桩-土体系相互作用弱化研究较少.

本文通过开展变水头渗透试验以及温控三轴试验对冻融条件下高原山区土体力学参数劣化规律进行研究，构建了基于土体结构劣化的水-热-力耦合模型并对其可靠性进行了验证.在此基础上，以实际工程为依托，对常年反复冻融循环条件下桩周土体以及桩基冻胀融沉特性进行了数值模拟与分析.

1 室内试验

1.1 试验土样

试验所用土样取自西藏日喀则萨迦县至秋洛村公路工程曲洛3 号中桥桩基钻孔，如图1 所示.依据《土工试验方法标准》（GB/T 50123—2019）［26］对土样进行室内土工试验，确定试样土样为粉质黏土，具有一定可塑性，其基本物理性质参数指标见表1.

图1 高原山区冻土Fig.1 Frozen soil in Plateau and mountainous areas

表1 原状土样基本物理性质Tab.1 Basic physical property of undisturbed soil

1.2 试验过程

1.2.1 变水头渗透试验

为探究冻融作用对土体渗透特性的影响，开展不同冻融循环条件下土样的变水头渗透试验.根据《土工试验方法标准》（GB/T 50123—2019）［26］制备渗透试验标准环刀土样.首先对试样进行冻融处理，依据日喀则当地工程地质资料［27］，设定冻结温度为-10 ℃，融化温度为20 ℃，冻结与融化时间均为12 h，以保证土样完全冻结与充分融化.整个冻融过程在TMS9018-250 型冻融循环试验箱（图2）中进行，最终完成制备0 次、1 次、2 次、4 次、6 次、8 次、10次冻融循环后的试样.

冻融试样制备完毕后，将其置于真空饱和器中进行抽真空饱和2 h，待无气泡产生后表明土样已达到饱和状态.将饱和后的试样放入变水头渗透试验装置中，试验装置如图3 所示.首先记录起始水头高度，每间隔5 min 观测记录水头高度变化，每组试样连续计测3次.每组试样重复试验5次后取平均值作为其渗透系数，冻融循环条件下的渗透试验方案如表2所示.

图3 变水头渗透试验装置Fig.3 Variable head permeability test device

表2 冻融循环条件渗透试验方案Tab.2 Penetration test scheme under freeze-thaw cycle conditions

1.2.2 温控三轴试验

为研究冻融作用对土样力学特性的影响，开展不同冻融循环条件下土样的温控三轴试验.首先对制备好的试样进行冻融循环处理，在冻融循环箱中完成1 次、2 次、4 次、6 次、8 次、10 次冻融循环后将试样取出，采用保鲜膜密封以防止冻融过程中水分的蒸发.

三轴试验采用GDS 液压动三轴测试系统（ELCTTS-2017）进行，如图4 所示，该设备主要由温控装置、剪切装置和控制系统三部分组成.其可施加最大轴力为10 kN，最小轴力为0.001 kN，所产生的最大法向位移为90 mm，最小位移为0.001 mm，围压最大可施加至70 MPa.经过不同冻融循环后的试样分别在不同有效围压下进行固结不排水（CU）三轴剪切试验，设定固结不排水剪切速率为0.07 mm/min，试验方案如表3所示.

图4 GDS液压动三轴测试系统（ELCTTS-2017）Fig.4 GDS hydraulic triaxial test system（ELCTTS-2017）

表3 冻融循环条件下温控三轴剪切试验方案Tab.3 Temperature controlled triaxial shear test scheme under freeze-thaw cycle conditions

1.3 结果与分析

1.3.1 冻融循环对土体渗透特性的影响

冻融循环作用对土体渗透特性具有重要影响［12，28］.试验土体渗透系数KT可通过式（1）获得.

式中：A为试样的横截面积，cm2；a为变水头管的横截面积，cm2；L为试样高度，cm；t1、t2分别为测读水头高度的起始、终止时间，s；h1、h2分别为起始水头和终止水头的高度，cm.

图5 表示土体渗透系数随冻融循环次数变化关系.由图5 可知，土体渗透系数随冻融循环次数的增加呈对数型上升趋势.未经过冻融循环（N=0）土样，初始渗透系数KT0=4.46×10-6cm/s；1 次冻融循环（N=1）后，土体渗透系数增大，达到6.36×10-6cm/s，较未冻融处理土样（N=0）增大了42.6%；6 次冻融循环（N=6）后，土体渗透系数达到1.05×10-5cm/s，较未冻融处理土样增大了135%.此后，当N＞6 时，土体渗透系数变化趋于稳定，最终经历10 次冻融循环后的土体（N=10），渗透系数达到1.2×10-5cm/s，较未冻融处理土样（N=0）增大了169%.这主要是由于冻融循环作用下土颗粒作为土体骨架发生了不可逆的损伤破坏，导致土颗粒间的孔隙增大.

图5 渗透系数与冻融循环次数变化关系曲线Fig.5 Variation curve of permeability coefficient and freeze-thaw cycle

不同冻融循环条件下土体孔隙比可通过表征软黏土孔隙比与渗透系数的Taylor公式［29］求得：

式中：e和e0分别为土体的当前孔隙比和初始孔隙比；k和k0分别为土体的当前渗透系数和初始渗透系数.其中土体转换系数Ck=0.5e0，土样的初始孔隙比e0=0.25.

图6 表示土体孔隙比与不同冻融循环次数之间的关系曲线.分析可知，土体孔隙比随冻融循环次数的增加呈对数型上升趋势.N=0 时，土体的初始孔隙比e0为0.249；1 次冻融循环（N=1）后，土体的孔隙比增大至0.254，较未冻融处理土样增大了2%；经历6次冻融循环（N=6）后，土体的孔隙比增大至0.261，较未冻融处理土样（N=0）增大了5%；当试样最终经历10 次冻融循环过程后（N=10），土体的孔隙比增大至0.267，较初始土样（N=0）增大了7%.产生该现象的主要原因为，土体冻结过程中，孔隙水由液态相变成冰，其体积增大，这一过程中土体内部由于相变产生的应力作用使土体内的孔隙增大；而在土体融化的过程中，土体内的冰晶重新相变成液态，但此时土体骨架来不及收缩或已发生不可逆转的损伤，因而总体来看土体内的孔隙是在不断增大的.

图6 孔隙比与冻融循环次数关系曲线Fig.6 Relationship curve between void ratio and freeze-thaw cycles

1.3.2 冻融循环对土体剪切特性的影响

图7 表示不同冻融循环条件下土体固结不排水（CU）三轴剪切试验过程中应力-应变关系曲线.由图7 可知，土体的抗剪强度随轴向应变的增大逐渐增大，均呈应变硬化特征.相同冻融条件下，围压对土体抗剪强度特征影响突出.以6 次冻融循环后的土体应力-应变曲线为例（图7（c）），取轴向应变ε=15%时所对应的剪应力为土体的破坏强度.当σ3=100 kPa 时，试样破坏强度为152 kPa；当σ3=200 kPa时，试样破坏强度为221 kPa，增大了45.4%；当σ3=300 kPa 时，破坏强度则为335 kPa，较σ3=100 kPa时，其破坏强度增大了120.4%.值得注意的是，试样经历6 次冻融循环后，其轴向应变发展出现了明显的超前特征.在剪应力还未发挥剪切作用时，试样就已经产生了变形，且随着冻融循环次数的增加，该特征更显著.

图7 不同冻融循环次数下三轴剪切应力-应变关系曲线Fig.7 Triaxial shear stress-strain curve under different freeze-thaw cycles

为更直观分析冻融循环次数对土体破坏强度的影响，绘制试样破坏强度与冻融循环次数之间的关系曲线，如图8 所示.土体破坏强度随冻融循环次数增加呈现不断减小的趋势.当σ3=100 kPa 时，10 次冻融循环后的试样（N=10）相较于未经过冻融循环试样（N=0）的破坏强度降低了40.2%；当σ3=200 kPa时，10 次冻融循环后的试样（N=10）相比于未经过冻融循环试样（N=0）的破坏强度降低了32.1%；而当σ3=300 kPa 时，这一数值则降低了25.7%.寒区土体在冻融循环的过程中其强度不断降低，土体破坏强度受冻融循环影响不断减小，这表明土体结构在冻融循环过程中产生了明显的劣化.

图8 土体破坏强度随冻融循环次数的变化Fig.8 Variation of soil failure strength with the number of freeze-thaw cycles

根据摩尔-库仑理论，可计算得到不同冻融循环次数下土体破坏强度对应的黏聚力c、内摩擦角φ，汇总结果如图9、图10所示.

图9 黏聚力与冻融循环次数变化关系曲线Fig.9 Relationship curve between cohesion and freeze-thaw cycle times

图10 内摩擦角与冻融循环次数变化关系曲线Fig.10 Relationship curve of internal friction angle and freeze-thaw cycles

图9 表示土体黏聚力随冻融循环次数变化关系曲线.可见，土体黏聚力随冻融循环次数的增加呈指数型下降.当N=0 时，土体的初始黏聚力c0=28.96 kPa；经历了1 次冻融循环（N=1）土体的黏聚力降至25.79 kPa，较初始试样（N=0）降低了11%；经历了6次冻融循环后（N=6）土体黏聚力降至14.45 kPa，较初始试样（N=0）降低了50.1%；而当经历了10次冻融循环后（N=10），土体的黏聚力最终降至11.05 kPa，较未冻融处理试样降低了61.8%.结果表明，第1 次冻融循环后试样黏聚力下降幅度最为明显，随着冻融循环次数的增加，前6 次冻融循环造成的累积效应对试样黏聚力的影响较大；后续随冻融次数的继续增加，对土体黏聚力的影响逐渐减弱.其原因主要是孔隙水冻结成冰导致土体体积膨胀，密度减小，土颗粒间的联结方式发生变化，导致土体黏聚力不断降低.随着冻融循环过程的继续发生，有部分土体持续发生着结构变化，但一部分原生土体结构已发生不可逆转的破坏，在此过程中表现为土体黏聚力的变化速率不断减小.

图10 表示内摩擦角与冻融循环次数变化关系曲线.分析可知，土体的内摩擦角随着冻融循环次数的增加呈对数型缓慢增加.N=0时，土体的初始内摩擦角φ0=16.37°；1 次冻融循环过程后（N=1），土体的内摩擦角为17.31 °，较未冻融处理试样（N=0）增大了5.7%；经历了6 次冻融循环后（N=6），土体的内摩擦角为19.07°，较未冻融处理试样增大了16.5%；10次冻融循环后（N=10），土体的内摩擦角为20.67 °，较未冻融处理试样（N=0）增大了26.3%.土体在冻融循环过程中，内摩擦角逐渐增大但其增长速率却在逐渐降低.这是因为孔隙水凝结成冰，孔隙体积增大，土颗粒间的有效接触面积减小，内摩擦角增大.随着冻融循环过程的继续，前一次融化过程中土颗粒深层孔隙中的冰晶体未完全融化，再一次冻结过程中孔隙水结成冰的相变过程较前一次有所削弱，最终表现为内摩擦角的增大速率逐渐降低.

通过对比分析冻融循环过程中黏聚力和内摩擦角的变化幅度，10 次冻融循环后土体黏聚力降低了61.8%，而土体内摩擦角仅增大了26.3%.表明冻融循环对内摩擦角的影响远远小于其对土体黏聚力的影响.

2 基于土体劣化桩-土体系数值模拟

2.1 基于土体结构劣化水-热-力耦合数学模型

2.1.1 土体渗透特性劣化特征

由1.3.1 小节分析可知，在冻融循环作用下，寒区土体的渗透系数以及孔隙比呈对数型增大.以冻融循环次数N为变量，构建渗透系数KT、孔隙比e值的劣化模型，对渗透系数KT值和孔隙比e值的劣化系数、Xe进行定义，如式（3）（4）：

式中：、Xe分别表示渗透系数、孔隙比的劣化系数；KT0、e0分别表示未经过冻融循环处理试样的渗透特性参数值和孔隙比参数值；、en分别表示n次冻融循环后试样渗透特性参数值和孔隙比参数值.

对不同冻融循环次数条件下渗透系数和孔隙比的劣化系数运用广义最小二乘法求得劣化系数的拟合曲线方程，拟合结果如图11所示.

图11 土体渗透特性参数劣化系数拟合曲线Fig.11 Deterioration coefficient curve of soil permeability parameters

由此可推导得到N次冻融循环后寒区土体的渗透系数参数、en：

2.1.2 土体剪切特性劣化特征

由1.3.2 小节的试验结果分析可得，土体的黏聚力c随冻融循环呈指数型降低，土体的内摩擦角φ随冻融循环呈对数型增加.以冻融循环次数N为变量，构建黏聚力c、内摩擦角φ的劣化模型，定义c、φ的劣化系数Xc、Xφ分别为：

式中：Xc、Xφ分别表示土体黏聚力、内摩擦角的劣化系数；c0、φ0分别表示未经过冻融循环处理试样的黏聚力以及内摩擦角参数值；cn、φn分别表示经过n次冻融循环后试样黏聚力以及内摩擦角参数值.

将不同冻融循环次数条件下土体黏聚力和内摩擦角的劣化系数进行计算汇总，并运用广义最小二乘法进行劣化系数曲线方程的拟合，拟合结果如图12所示.

图12 土体剪切特性劣化系数拟合曲线Fig.12 Deterioration coefficient curve of soil shear properties

由此可推导出N次冻融循环后土体的黏聚力cn以及内摩擦角φn表达式分别为：

2.1.3 土体水-热-力耦合控制方程修正

土体水-热-力耦合控制方程的修正，是在归纳总结国内外学者们的研究方法与理论基础上［30］，充分考虑土体劣化在桩服役期间的影响，对现有的控制方程进行的改进.

1）水分场控制方程

根据质量守恒定律建立的水分场控制方程：

式中：θl为未冻水体积分数；θi为冰体积分数；ρ1为水的密度，g/cm3；ρi为冰的密度，g/cm3；v为水流通量；ψm为基质势；z为重力势；KT（θl）为土体的渗透系数，cm/s.

由2.1.2 节所得结论，土体的渗透系数同时也是冻融循环次数N的函数，将式（5）得到的土体渗透系数参数KTn，代入水分场控制方程进行修正得：

2）温度场控制方程

根据能量守恒定律建立温度场控制方程：

式中：C为土体热容，kJ/（m3·K）；L为相变潜热，kJ/kg；T为温度，℃；λ为导热系数，W/（m·K）.

在冻融循环过程中，温度场与水分场之间相互影响，具有强耦合关系.温度场控制方程与水分场控制方程耦合过程中存在θl、θi以及T三个未知函数，因此需要一个联系方程实现水热两场耦合并进行修正.通过引入白青波等［31］定义的固液比B（T）作为两场之间的联系方程：

式中：Tf为土体的冻结温度，℃；b为固液比系数，粉土取值0.47，黏土取值0.56.

3）应力场控制方程

由水-热两场耦合产生的体积应变作为土体内部应力场建立的控制方程如式（15）：

式中：θ0为初始含水量；Δθ为水分迁移量；θu为未冻水含量；e为初始孔隙比.

在冻融循环过程中，孔隙比e随着冻融循环次数N的增加不断变化，将式（6）代入式（15）后可得修正后的应力场控制方程：

2.1.4 模型验证

基于时伟等［32］研究冻融循环作用下膨胀土动力学特性的土柱试验进行数学模型验证.土柱直径为70 mm，高140 mm，如图13 所示.由于实际模型试验中对土柱的底部和四周采取了保温措施，因此在构建模型时，假定土柱与外界之间的热量交换仅从顶部发生.水分场边界条件设置为零通量，温度场边界条件设置6 次单向冻融循环试验，数值模型计算所需参数如表4所示.

图13 土柱单元模型Fig.13 Mesh generation of soil column element model

表4 土柱数值模型参数Tab.4 Parameter table of soil column numerical model

图14 为土柱顶部位移随时间的变化曲线.由图14 可知，随着冻融循环次数的增加，冻融作用对土柱的冻胀融沉特性的影响开始减弱，这与前文中冻融作用对土体的渗透特性与剪切特性参数值的影响规律一致，证明本文所提出的劣化模型的合理性.同时，考虑了土体结构劣化的模拟值更接近试验值，且模拟值与实测数据误差在5%以内.而未考虑结构损伤模型计算出的最终冻胀位移要比实测值小0.41 mm，误差约为28%.进一步表明了基于土体结构损伤的水-热-力耦合数学模型的可靠性.

图14 位移随时间变化曲线Fig.14 Displacement versus time curve

2.2 桩-土体系数值模拟

2.2.1 工程概况

以日喀则萨迦县至秋洛村公路工程曲洛3 号中桥桩基为实际工程背景，利用COMSOL 有限元数值软件中系数型偏微分方程模块进行二次开发，构建二维桩-土体系数值模型，对桩基服役50 年内桩-土体系位移场进行预测分析.曲洛3 号中桥桩长为20 m，桩直径为1.2 m，选用钻孔钢筋灌注桩，混凝土标号为C30，桩周土体以粉质黏土为主.

2.2.2 模型建立

基于2.1.3 节中建立的考虑土体结构劣化的水-热-力三场耦合数学模型，利用有限元软件COMSOL系数性偏微分方程进行二次开发，将温度场、水分场以及应力场控制方程转换为系数性偏微分方程形式，如式（17）所示.设置桩周土体材料，将土体剪切特性劣化参数输入.

图15 为二维桩-土体系数值模型，设置桩长为20 m，桩径1.2 m，通过试算设定土体的厚度为30 m，土体宽度为12 m（10倍桩径），使用自由三角形划分网格，桩-土界面使用自定义接触对，以固定节点定义接触.桩体及桩周土的相关热力学参数如表5、表6所示.

表6 桩周土热力学参数Tab.6 Thermodynamic parameters of soil around pile

图15 桩-土体系有限元数值模型Fig.15 Finite element numerical model of pile-soil system

表5 桩体热力学参数Tab.5 Thermodynamic parameters of pile

设置桩-土体系数值模型左右边界为绝热条件，无热传导.上边界设置为当地大气温度，预测50 年后平均气温升高1～2.6 ℃［33］，定义温度解析函数；同时，上边界采用狄式边界条件，设置水分场边界为无通量，上边界为自由边界；左边界为对称边界，右边界为辊支撑，限制左右位移，底部为固定约束.

2.2.3 预测结果分析

图16 为桩周土最大位移随时间变化曲线.分析可知，随着桩基服役年限增加，土体的最大冻胀位移不断减小，冻胀速率变慢.桩基服役前期，随着冻融循环的进行，桩周土体原生结构破坏严重，土体上部含水量增长较快.冻结条件下，液态水的相变导致土体产生冻胀；桩基服役后期，土体结构受冻融循环影响变小，上层含水量增量减小，冻胀位移变化趋于平缓.另外，从图16 中可看出，桩周土体的最大融沉量随桩基服役年限的增加逐渐增大，其主要原因为土体结构劣化导致冻融层含水量增加.与桩周土体冻胀位移变化规律类似，冻融循环后期冻融作用对土体融沉位移影响减小.整体来看，融沉位移要远大于冻胀位移，桩周土体呈下沉趋势.该下沉趋势对桩基产生一定的负摩阻力，同时桩-土界面冰膜的融化会进一步降低界面剪切强度.因此，土冻融循环会削弱桩基承载力，对长期服役桩基承载性能带来不利影响.

图16 桩周土体最大位移随时间变化Fig.16 Variation of soil displacement around pile with time

图17 为桩基位移随服役年限的变化曲线.由图17可知，随着桩基服役年限增加，桩基冻胀位移变化趋于平缓，冻胀速率变慢.一方面，主要是因为土体结构劣化受冻融作用影响逐年降低；另一方面，桩基融沉量随服役年限的增加逐渐增大.该现象同样是由于土体结构劣化受冻融作用影响随时间减小，位移变化速率降低.总体而言，桩基服役期间，桩基融沉位移大于其冻胀位移，桩身整体呈融沉特征.对于高原山区实际工程而言，要正确看待桩基融沉现象并采取一定措施以减小其融沉位移的产生.