0 引言
在寒冷地区,大量工程实践表明,修建于冻土层上的构筑物会极大地改变冻土层原有的热平衡,进而引起多年冻土层的剧烈融化沉降,这将进一步导致构筑物基础的不均匀沉降破坏。冻土融化沉降会对多年冻土地区的铁路、公路、民用建筑和工程建筑等造成严重的危害[1-4]。冻土融化沉降实质是随着土体的不断融化,融化区域内的孔隙水在外压力作用下持续排出,固结变形不断发展的复杂物理力学过程[5]。在过去的几十年里,为了更好地揭示融化固结机理,学者们开展了大量的研究工作。Morgenstern等[6]基于太沙基固结理论和半经验融化边界方程建立了一维融化固结理论,但是该理论仅局限于一维情况,不能用于边界条件复杂的实际工程问题。针对这一缺点,Sykes等[7]将融化固结理论扩展至三维问题的运用,但是上述理论均基于小变形假设和线弹性理论,在计算高含冰量多年冻土时会产生较大的预测误差。为了提高预测精度,Foriero等[8]以孔隙比为自变量推导得到了一维大变形融化固结理论,不过该理论仍然局限于一维情况。考虑到边界条件复杂的实际工程问题,Yao等[9]基于欧拉描述进一步将融化固结理论扩展至计算三维大变形融化固结问题,并结合室内试验结果证明了其有效性。
上述的研究工作可以很好地描述具有恒定热边界的融化沉降固结特性,即在融化过程中,融化沉降不断发展,固结度保持恒定。而对于具有季节性变化热边界的实际工程案例的融化沉降问题[10-11],恒定热边界的融化固结理论不再适用。以青藏公路为例,路基融化沉降通常发生在暖季,而在冷季时停止,这意味着融化沉降不是全年发生的,而是随着路基表面温度上升排水边界融化开放而季节性变化发生的。并且高含冰量土体在融化后本构关系具有强烈的非线性特征[12],且强度较低,在较低的应力水平条件下也会产生一定的塑性变形[13-14]。故使用线弹性本构关系会进一步降低土体沉降的计算精度。因此,研究高含冰量冻土路基在季节性排水边界下的融化固结特性对工程应用具有重要意义。
本文将对已有的三维非线性大变形融化固结理论进行修正,引入季节性温度边界条件及摩尔-库伦破坏准则,建立了能够考虑季节性温度边界条件影响的三维非线性塑性融化固结理论模型,在此基础上,采用FLAC3D软件对所建理论模型进行数值实现。以青藏公路某段高含冰量路基为例,分析了其在季节性温度边界条件下的融化固结规律,并结合实测数据验证了所建理论模型的合理性。最后,进一步分析了冻土路基在融化固结过程中应力和孔隙水压力的分布及发展规律。
1 控制方程与数值实现
1.1 控制方程
式中:σij为总应力张量;σij'j为σij沿j方向的梯度;gi为重力加速度(m/s2);ρ为密度(kg/m3);vi、vj为瞬时速度(m/s),i、j代表不同质点; 、 均为旋转率张量; 为总应力率;Dijkl为应力-应变关系张量;δij为克氏记号; 为孔隙水压力变化率;xi与xj分别为物质点在t时刻的位置; 与 分别为应变率张量和旋转率张量。
假设土颗粒及其孔隙中的流体不可压缩,则流体的质量守恒方程可以表示为
式中:qi为孔隙水相对土体骨架的流速(m/s);qi'i为流体沿i方向的流速梯度;qv为体积流源(s-1);εv=εii(i=1,2,3)为体积应变;qi可用达西定律描述为
式中:k为介质渗透系数(m/s);u为孔隙水压力(kPa);ρw为流体密度(kg/m3);gj为重力加速度(m/s2)。
式中K和G分别为排水条件下的体积模量(Pa)和剪切模量(Pa)。可以用杨氏模量E与泊松比v表示为
而杨氏模量E可以进一步表示为
式中:Es为压缩模量(MPa);P'为有效应力(kPa);εv为体积应变。其中P'与εv可以用孔隙比e表示为:
式中:e0为初始孔隙比;β为通过K0固结试验获取的 e-lnP'曲线的斜率[15]。把式(11)—式(13)代入式(10)中,压缩模量Es可进一步写为
将式(14)代入式(7)中,用压缩模量与孔隙比的关系描述土骨架在融化固结过程中的非线性力学行为。
对于土体进入塑性阶段的沉降变形使用摩尔-库伦准则进行描述,具体如下:
式中:f为屈服函数;σ1和σ3分别为第一主应力和第三主应力;φ为摩擦角(°);c1为黏聚力(Pa);Nφ为与土体摩擦角φ相关的无量纲参数。
式中:T为温度(℃);hv为流体热源(W/m3);hi和hi'i分别为沿i方向的热流密度(W/m2)和热流密度梯度;T'i为沿i方向的温度梯度(℃/m);c与λ分别为考虑冰水相变的比热容(J/(kg·℃))和热传导系数(W/(m·℃)),与温度有关,可以表示为:
式中:cu、cf分别为融土和冻土的比热容(J/(kg·℃));λu和λf分别为融土和冻土的热传导系数(W/(m·℃));L为水的相变潜热(J/kg);W为含水量(%);Wi为含冰量(%);Tp和Tb分别为冰水相变的上下温度界限(℃)。
1.2 数值实现
冻土融化固结的计算问题实际上是热-水-力三场的耦合计算,是一种典型的循环计算。其包括2个阶段,即热学计算阶段与固结计算阶段。在每个循环过程中,首先在整个计算区域进行一个热学时间步长(dtth)的热学计算,根据计算结果确定融化区域(区域温度>0 ℃),并判断排水边界是否开放(边界温度>0 ℃)。然后,在已经融化的区域内进行n步流体时间步长(dtfl)的固结计算。其中,热学时间步长与流体时间步长必须满足关系dtth=n·dtfl,以保证每个循环结束时流体计算和热学计算的计算时间一致。在季节性温度边界的作用下,先前的融化区域会再次冻结,设置其渗透系数为0,此时固结过程不再发生;而当该区域再次融化时,设置其渗透系数为对应土层实际值,此时固结过程相应进行。具体计算流程如图1所示。
2 案例研究
2.1 现场说明
修建于20世纪50年代运营至今的青藏公路期间经历了多次大规模改建[16]。由于早期路基地温数据缺失,加之公路沿线地表温度条件随工程改建发生了显著变化,导致采用数值方法对路基融化固结过程进行实时模拟面临极大困难。因此,本文将路基近年来天然孔所监测的温度作为初始温度条件,将路基中心、路肩和天然地表处距离地表最近的监测点温度作为边界条件(距离地表0.3 m),模拟选定的多年冻土路基50 a来的融化固结过程。根据路基的工程地质条件可以确定各层的热学和力学参数,如表1所示[5,16 -18]。文献[12]—文献[14]的研究表明,对于高含冰量冻土的融化固结问题,线弹性应力-应变关系会造成较大的误差,而对于低含冰量冻土,其误差较小。因此,为了提高计算效率,本文仅对黏土层和亚黏土层2层高含冰量土层进行非线性塑性应力-应变关系分析,而对路基填土层和强风化泥岩层两层低含冰量土层进行线弹性应力-应变关系分析。
表1 数值计算的热学与力学参数Table 1 Thermal and mechanical parameters for numerical calculation |
| 岩性 | 岩层 范围/ m | 密度/ (kg·m-3) | 含水 量/ % | 比热容/ (J/(kg·℃) -1) | 热传导系数/ (W/(m·℃) -1) | 压缩 模量/ MPa | 泊松 比 | 渗透 系数/ (m·s-1) | 黏聚 力/Pa | 内摩 擦角/ (°) | 初始 孔隙 比 | β | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 融土 | 冻土 | 融土 | 冻土 | |||||||||||
| 路基 填土 | 0~3 | 2 230 | 10 | 2 143 | 1 627 | 1.31 | 1.71 | 1.65 | 0.25 | 1.00×10-7 | 3.50×104 | 22 | — | — |
| 黏土 | 3~6.5 | 1 810 | 50 | 1 894 | 1 283 | 1.18 | 1.86 | — | 0.30 | 1.00×10-9 | 3.30×104 | 25 | 1.35 | 0.080 |
| 亚黏土 | 6.5~9.5 | 2 150 | 40 | 1 397 | 1 138 | 1.21 | 1.52 | — | 0.25 | 1.00×10-9 | 3.00×104 | 24 | 1.08 | 0.052 |
| 强风化 泥岩 | 9.5~30 | 2 160 | 15 | 2 191 | 1 907 | 0.91 | 0.96 | 5.00 | 0.25 | 5.00×10-9 | 7.00×105 | 48 | — | — |
2.2 边界条件
式中:T0为年均温度(℃);θ为气温变暖速率(0.02 ℃/a)[19];l为温度振幅(℃);t为时间(s);nπ为初始相位,相关参数取值见表2[17]。路基的初始温度场根据图4给出的路基天然孔监测温度确定。从夏季温度最高值开始计算,即路基中心和路肩的温度初始相位为π/2。路基中心孔和路肩孔的温度峰值出现于同一时间,天然地表温度峰值出现时间较路基中心孔和路肩孔的温度峰值滞后约一个月,故天然地表的初始相位为π/3。路基填土的初始温度取为6 ℃。本案例仅考虑路基在自重作用下的固结变形。路基模型的上边界设置为自由位移边界,左、右边界不沿水平方向移动。土体冻结时沿竖直方向的位移固定不变,在融化时设置为自由边界。沥青路面设置为不透水边界,路基边坡和天然地表在融化时设置为自由排水边界。
3 结果与分析
3.1 融化深度和变形
通过对比图5中路基融化深度的计算结果与实测结果可以发现2004年的实测结果与第25年的计算结果基本一致。故本文将融化沉降第25年至第27年的计算结果与2004—2006年的实测结果进行比较分析。
图6给出了路基中心考虑塑性变形与不考虑塑性变形2种融化沉降计算结果与实测结果的对比。由图6可知,路基中心融化沉降的2种计算结果和实测结果随时间的发展表现为暖季融化沉降冷季冻结停滞的周期性变化规律,即在每年暖季到来路基排水边界(路基边坡和天然地表)融化开放时,路基的排水通道随之开放,土体中的孔隙水逐渐排出,土骨架发生融化沉降,即固结过程。固结过程持续进行直至冷季来临。冷季来临之后,由于地表土体冻结(路基排水边界关闭),排水路径被冻结封堵,土体固结过程也随之停止,路基沉降的发展也停滞。进一步对比考虑塑性变形与不考虑塑性变形两种计算结果和实测结果可以看出,相较于不考虑塑性变形的计算结果,考虑塑性变形的计算结果更接近于实测结果,并且由于考虑塑性变形所产生的影响随着时间的推移不断积累,故长期计算高含冰量冻土工程融化固结变形时,塑性变形所产生的影响不可忽略。尽管考虑塑性变形的计算结果更接近于实测结果,但与实测结果仍有较大差距。
3.2 应力场
为了研究高含冰量冻土路基融化固结过程中的应力传递及分布规律,分析了融化沉降后区域的垂直有效应力场。图7给出了10 a与30 a的冻土路基内垂直有效应力分布的等值线。研究表明[5,12],冻土路基表层路基填土的融化沉降过程通常在运营初期(<5 a)内完成,后续冻土路基的融化沉降主要是更深的永久冻土层或新融化层的融化固结。此外,根据10、20、30 a的冻土路基内部孔隙水压力分布的计算结果表明(图8),冻土路基内部的孔隙水向下并由路基边坡及天然地表向外排出,即在运营初期之后,路基填土中已经没有孔隙水压力。相应地,本文计算的路基填土中有效应力近似等于总应力,图7给出的路基内部有效应力的分布可用于描述路基的应力状态。因此,在图7中,采用计算得到的垂直有效应力来表示冻土路基融化后区域的应力状态。由图7可知,计算得到的10 a与30 a的垂直有效应力等值线的发展规律相似,即不同时间的垂直有效应力等值线均呈锅形分布,这主要是受路基土层融化区域融化深度分布的影响(图5)。
进一步对比10 a与30 a的垂直有效应力等值线可知,随着时间的发展,垂直有效应力等值线的分布范围不断扩大,相同等值线所围成的面积也在逐渐增大(如120 kPa垂直有效应力等值线所围成的面积)。这主要是由于随着时间的发展,融化深度在逐渐增大,相应地,由于路基土体自重的影响,垂直有效应力也在逐渐增大。
3.3 孔隙水压力
冻土路基的固结过程主要受孔隙水压力分布的影响[15]。图8给出了10、20、30 a的冻土路基内孔隙水流向以及孔隙水压力分布的等值线,其中,箭头表示冻土路基内融化区域孔隙水流动方向。如图8所示,在不同的时间,由于路基自重和边界条件(沥青路面不透水)的作用,冻土路基融化区域内的多余孔隙水由路堤边坡和天然地表处排出。根据冻土路基内孔隙水压力分布的等值线可以看出,大部分孔隙水压力主要作用于最大融化深度的附近区域,而在路基浅层区域,大部分孔隙水压力在10、20、30 a的时间已经完全消散,也就是说,路基浅层区域融化后的孔隙水在运营初期已经由路堤边坡和天然地表处排出。在之后的长期运营过程内,冻土路基融化沉降的发展主要是由于新融化的孔隙水在最大融化深度处或永久冻土区域上限处的消散。这是多年冻土路基融化沉降规律与未冻土层路基之间的主要差异之一。
4 结论
(1)本文建立了能够考虑季节性温度边界条件影响的三维非线性塑性融化固结理论模型,基于该理论模拟研究了青藏公路某段高含冰量冻土路基的融化固结问题,通过将计算结果与实测数据进行对比分析表明,改进的理论模型能够更合理地描述季节性温度边界条件下高含冰量冻土路基的融化固结特性。
(2)冻土路基的沉降随着地表温度的季节性变化而呈现出暖季融化沉降而冷季冻结停滞的周期性变化规律,这是季节性温度边界条件下冻土路基融化固结规律的最显著特征。
(3)高含冰量冻土路基融化跑龙套结过程中塑性变形随时间持续累积,因此在高含冰量冻土路基长期变形预测中,塑性变形所产生的影响不可忽略。
(4)冻土路基浅层融化区域内的孔隙水在运营初期已经消散,而在之后长时间的运营过程中,冻土路基融化沉降的持续发展主要是由于融化锋面处新融化的孔隙水的消散,这是多年冻土路基融化沉降规律与未冻土层路基的主要差异之一。