寒区非圆形隧道冻胀力的解析解1)

2020-02-23李岩松陈寿根

力学学报 2020年1期

李岩松陈寿根

∗(中国五冶集团有限公司,成都 610063)

†(西南交通大学交通隧道工程教育部重点实验室,成都 610031)

引言

在高寒高海拔地区修建隧道,隧址区原有的热力平衡被打破,隧道的通风对流使冷空气在隧道内部的流动,造成了隧道内部温度的下降,从而创造了隧道围岩形成季节冻土或多年冻土的条件.当围岩温度过低围岩中的裂隙、孔隙水冻结时会发生体积膨胀,在这种体积膨胀受到隧道衬砌和未冻岩体的约束时,冻结围岩便会对隧道衬砌产生冻胀力.20 世纪70∼80 年代,日本对东北寒冷地区隧道的调查分析结果表明,冻胀力是造成寒区隧道变形破坏主要因素.由于我国经济的高速发展,高速铁路、公路建设迎来了高速发展时期,不可避免地需要在高寒高海拔地区修建隧道.隧道常常发生冻害,对隧道的维修养护造成了极大困难,更是对隧道结构和运营安全构成了严重威胁.为保证隧道结构和运营安全,降低维修养护成本,各国学者从不同方向、不同领域对寒区隧道冻胀力进行了研究.

Lai 等[1-4]应用无量纲微扰法、Galerkin 法、黏弹性理论以及CT 技术分别对寒区隧道进行多项研究并取得了显著成果.为寒区隧道施工、运营、维护提供了有力的参考依据;Gao 等[5]应用经典弹塑性理论得到了圆形寒区隧道围岩应力和塑性区,发现冻胀力会导致围岩塑性区显著膨胀并在屈服后出现塑性硬化现象;Feng 等[6]应用弹塑性理论得到了圆形寒区隧道围岩应力及位移的弹塑性解析解,并且通过参数分析,得出相关参数对冻结圆周内、外边缘塑性半径和围岩应力的影响结果;夏才初等[7]考虑施工过程中机械产生的热量对多年冻土的融化影响,分析了多年冻土融化范围,建立并求解融化作用下多年冻土隧道围岩弹塑性模型,确定了相应控制融化作用的关键因素;张常光等[8]在弹性模量变化、中间主应力效应、围岩应变软化和剪胀条件下,基于统一强度理论推导了深埋圆形隧道塑性位移新解;谭贤君[9]、杨罡[10]、黄继辉等[11-12]、苏林军[13]应用有限元软件对冻土隧道冻胀力进行模拟,得到了一系列令人满意的结果;吴剑[14]通过模型试验结果与数值模拟结果对比,详细研究冻土隧道冻胀力,得出了不同冻胀条件下隧道衬砌变形特征曲线以及冻胀力分布规律;根据Takashi 方程Hao 等[15]得到了圆柱形冻土试块的冻胀力以及温度分布规律,论证了该方法在三维空间中的适用性;Tan 等[16]应用数值模拟建立模型,计算了寒区隧道在多孔介质中冻融循环的温度扩散问题,通过数值模拟得到隧道衬砌表面的建筑保温材料会极大地改变衬砌和围岩温度的结论.

在以上寒区隧道研究中所涉及应用弹性、弹塑性理论计算寒区隧道冻胀力解析解时均采用圆形隧道模型,并未采用与实际隧道断面相似的隧道模型进行计算.由于计算模型的差异,这些研究成果与实际结果存在一定的差异,不能直接应用于非圆形寒区隧道问题中.

当考虑支护为圆形的条件时,寒区隧道冻胀力和位移的弹性解析解通过连续性条件易于给出.而当隧道断面形状为马蹄形、直墙拱形等非圆形状并且考虑隧道支护时,由于考虑了隧道支护,使问题由单连通域问题变为多连通域问题,不能直接应用经典的复变函数理论求解非圆形隧道应力和位移.所以,还鲜有学者给出考虑衬砌因素的非圆形寒区隧道冻胀力和冻胀位移的弹性解析解.综上所述,为得到考虑衬砌的非圆形寒区隧道冻胀力弹性解析解,应用经典复变函数理论幂级数解法,结合连续性条件,在文献[17-18]的基础上,推导了考虑衬砌的非圆形寒区隧道冻胀力的弹性解析解.为下一步非圆形寒区隧道冻胀力的弹塑性分析奠定了基础.

1 问题描述

根据弹性理论,寒区隧道在深埋条件下,由于隧道开挖尺寸与埋深相比很小,隧道开挖问题可视为在无限平面内的孔洞问题,不受边界效应影响,地应力变化梯度可忽略不计.当考虑冻胀圈和支护结构时,无限平面内的孔洞问题由单连通域孔洞问题变为多连通域问题,无限平面被分割为3 个区域.图1 给出了无限平面内考虑冻胀圈的隧道孔洞模型.其中在Z平面中衬砌、冻胀圈和未冻围岩分别表示为S1,S2和S3区域,所对应的边界L1为衬砌与大气接触面,它为自由边界,受到的面力为零;L2为衬砌与冻胀圈接触面,它满足连续边界条件;L3为冻胀圈与未冻围岩接触面,它满足连续边界条件;围岩的远场应力由隧道埋深边界条件确定.为了进行计算,必须将Z平面上的非圆形衬砌、冻胀圈结构通过保角变换映射到ζ 平面上,使得在Z平面上的非圆形结构变换为ζ平面上的圆环结构.这时Z平面上的区域S3映射为ζ 平面行的区域O3;区域S2映射为ζ 平面行的区域O2;衬砌区域S1映射为圆环区域O1;边界L1,L2及L3分别映射为圆环R1,R2与R3.其中最一般的映射函数可写成洛朗级数的形式,即

图1 Z 平面保角变换为ζ 平面示意图Fig.1 Conformal mapping of tunnel in z-plane in to two concentric circles in ζ-plane

其中,R是一个实常数,它的大小反映了隧道孔洞的大小;Ck一般为复常数,其中k=1,2,3,···,Ck反映了隧道断面形状,在一般情况下Ck取前几项时,映射函数的精度已经可以满足工程要求[19].ζ 为Z平面上的点映射到ζ 平面上的正交曲线坐标,通常可由角度θ 与半径ρ 表示.

2 复变函数解法公式推导及求解

根据Kargar 等[20-22]和陈子荫[23]提出的方法,ζ面上未冻围岩O3区域内的应力函数可表示为

式中,hk和mk为所求复常数,它们均由边界条件确定;φ0(ζ)和ψ0(ζ)满足φ0(∞)=0 和ψ0(∞)=0;Lu等[21-22]提出Γ 和Γ′可由无穷远处应力边界条件定出,即

在方程(6)和(7)中K为侧压力系数,γ 和H分别表征围岩重度与隧道埋深.

ζ 平面上冻胀圈O2区域的应力函数φ2(ζ)和ψ2(ζ)可表示为

其中,ak,bk,ck和dk为所求复常数,它们均由边界条件确定.

ζ 平面上衬砌O1区域应力函数φ3(ζ)和ψ3(ζ)可表示为

其中,Ak,Bk,Ck和Dk为所求复常数,它们均由边界条件确定.

在复变函数解中[24-26]应力均可以用两个应力函数φ 和ψ 表达,衬砌、冻胀圈和未冻围岩上应力均可表示为

其中,σρ,σθ和τρθ分别为正交曲线坐标下的法向应力、环向应力和剪应力.

假设寒区隧道衬砌S1和未冻围岩S3区域体积不受围岩温度变化影响,相应ζ 平面上映射区域O1和O3位移满足复变函数[24-26]应力函数表达,衬砌和未冻围岩上位移可表示为

其中,uρ和uθ分别表示法向位移与环向位移,G是切变模量,κ 由泊松比确定,κ=3-4υ 和κ=(3-υ)/(1+υ)分别为平面应变与平面应力状态下κ的取值,υ 为泊松比,因隧道纵向长度与断面最大直径相比很大,故问题可视为平面应变问题,κ 取3−4υ.

由于冻胀圈S2区域内体积受温度变化影响,故相应在ζ 平面上O2区域内位移无法满足复变函数[24-26]应力函数表达,需要重新推导(推导过程详见附录1),经推导冻胀圈上位移可表示为

其中ε 为冻胀圈线膨胀系数,可由冻胀区单位冻胀体积变化定出,即

其中∆V是冻胀区单位冻胀体积变化.

应力函数φ1(ζ),ψ1(ζ),φ2(ζ)和ψ2(ζ),φ3(ζ)和ψ3(ζ)必须分别在R2和R3上满足连续边界条件,在R1上满足应力边界条件.在衬砌与大气接触面上力的边界条件可表示为

在计算非圆形隧道围岩应力时,为能更好地反应实际隧洞变形位移,Lu 等[20-22]在连续性位移条件中引入了开挖卸荷因子.为了能更直观反映冻胀力对非圆形隧道衬砌的影响,在位移连续边界条件中并未考虑开挖卸荷因子,如需考虑开挖卸荷效应,只需在位移边界条件中添加符合实际工况的开挖卸荷因子即可.根据推导的位移公式,在冻胀圈与衬砌接触面上位移连续边界条件可表示为

在冻胀圈与衬砌接触面上力的连续边界条件可表示为

在未冻围岩与冻胀圈接触面上,根据推导的考虑冻胀作用的位移公式(附录1),位移连续边界条件可表示为

在未冻围岩与冻胀圈接触面上力的连续边界条件可表示为

其中,t1,t2和t3分别表示ζ 平面上R1,R2和R3边界上的点.

方程(2)和(3)中的表达式Γw(ζ)和Γ′w(ζ)表示围岩的初始地应力和初始位移的分量,所以在方程(19)中计算由于开挖引起的位移时应将其去除[27];方程(17)∼(20)在计算过程中假设接触面之间无滑动.

应用于李岩松等[17-18]幂级数解法,方程(16)可分别展开为正幂次方程组(21)和负幂次方程组(22)

其中,Lv可由方程(1)中R和Ck计算得出,ζ 为Z平面上的点映射到ζ 平面上的正交曲线坐标,可表示为ζ=ρσ,这里σ=exp(iθ);ζ 的共轭可表示为为便于计算,方程(1)可写为

方程(17)可分别展开为正幂次方程组(24)和负幂次方程组(25);方程(18)可分别展开为正幂次方程组(26)和负幂次方程组(27);方程(19)可分别展开为正幂次方程组(28)和负幂次方程组(29),方程(24)∼(29)表达式如下

最后将方程(20)分别展开为正幂次方程组(30)和负幂次方程组(31).应当注意的是:由于边界条件不同,方程(21)和(22)中的极半径ρ 与方程(24)∼(27)以及方程(28)∼(31)中的极半径ρ 在数值上分别等于r1,1 和r2.

结合行列式(21)、式(22)以及行列式(24)∼式(31),可以得到无穷线性方程组[17-18],求解无穷线性方程组得到应力函数(方程(2)和(3),方程(8)和(9),方程(10)和(11))中各项系数.由得到的应力函数计算方程(12)和(13)得到非圆形寒区隧道应力解析解;由方程(14)计算衬砌和未冻围岩位移解析解;最后由方程(15)计算冻胀圈位移解析解.

3 工程案例应用分析

以新建汶川至马尔康公路鹧鸪山特长隧道洞口段为研究对象,利用上节解析方法,对非圆形寒区隧道应力和变形进行求解,并用Flac 有限差分软件进行进行校核,验证解析解的准确性.

根据现场大气、围岩实测温度、围岩及结构热力学试验,估算围岩冻结深度.可知当原岩温度为10◦C时,冻结深度约为2 m[28].据此确定围岩冻胀圈和未冻围岩.

结合鹧鸪山隧道地质报告、室内岩石试验和点荷载试验数据[28].围岩、结构材料物理力学参数见表1.

表1 地层及材料参数Table 1 Main physical parameter for tunnel calculation

图2(a)为鹧鸪山隧道K188+480 m 处断面计算简图.图2(b)为Flac 有限差分软件计算模型图,共划分网格3916 个,用于模拟深埋非圆形寒区隧道开挖后衬砌、冻胀圈和围岩的应力及变形状态.隧道采用全断面开挖,开挖后立即支护;假定模型为弹性模型;衬砌与冻胀圈,冻胀圈与围岩之间只受法向力;隧道冻胀系数采用等效热胀系数进行模拟;假设冻胀圈温度恒定,无温差;隧道冻胀力和冻胀变形通过控制热力模型和结构模型交替开关进行计算;衬砌与冻胀圈,冻胀圈与围岩之间无滑动.根据式(15)中对线胀系数的说明、水冻结为冰体积增大约9%及冻胀圈含水量计算得冻胀圈线膨胀系数ε≈0.002 5.

根据吕爱钟等[19]提出的最优化法求解映射函数的思路,编写相应MATLAB 计算程序,通过计算,得到映射函数

在工程计算中方程(1)取前几项就可满足工程精度要求[19],所以为提高计算效率,式(37)取映射函数(方程(1))的前四项;式(37)中ζ=ρσ,其中当ρ=1.225 时,表示映射函数所对应未冻围岩与冻胀圈接触面的边界条件;当ρ=1 时,表示映射函数所对应冻胀圈与衬砌接触面的边界条件;当ρ=0.863 3时,表示映射函数所对应衬砌与大气接触面的边界条件.

图2 (a)隧道计算简图(单位:m);(b)隧道计算模型Fig.2 (a)Tunnel calculation diagram(Unit:m)and(b)Finite mesh calculation mode

根据边界条件(方程(16)∼(20))编写MATLAB计算程序,求得寒区非圆形隧道冻胀应力及冻胀位移解析解,并将其与Flac 有限差分软件所求的数值解进行比较,比较结果如下所示.

由于Flac 有限差分软件计算应力为隧道衬砌及冻胀圈的竖向和水平向应力,这与解析解所求法向、环向应力坐标系统不统一,无法进行比较,故将Flac有限差分软件所求应力结果进行坐标转换后再进行对比.

图3 是衬砌内边界环向应力,由图中可看到,数值解与解析解吻合得较好;由于冻胀作用的影响衬砌内边界环向应力显著增大,在拱顶(α=0◦)到拱肩(α=40◦)这一扇形区域衬砌所受的环向应力明显增大;在拱脚(α=120◦)处,环向应力集中明显,拱脚处环向应力较未冻结时增大约25.9%;在拱底(α=180◦)附近30◦的扇形区域范围内衬砌所受环向应力增长显著.

图3 衬砌内边界环向应力Fig.3 Circumferential stress along the inner lining periphery

图4 为衬砌与冻胀圈接触面上法向应力,由图4可以看出:衬砌法向应力较未受冻胀力时有明显波动.由于冻胀作用,衬砌在α=60◦附近法向应力较未受冻胀力时增大约107.5%;在α=110◦附近冻胀后衬砌与未冻胀衬砌相比法向应力减小约82.4%;在α=140◦附近受冻胀力影响法向应力较未冻胀法向应力增大约21.4%;在衬砌与冻胀圈接触面上,数值模拟所求法向应力与解析法所求法向应力吻合得较好,只是在部分路径上有一定的偏差.

图4 衬砌与冻胀圈接触面上法向应力Fig.4 Normal stress along frost heaving ring-lining interface

由图5 可见,冻胀作用对衬砌影响明显,尤其在衬砌拱顶、拱脚、拱底处.由于冻胀力影响,拱顶环向应力由2.01 MPa 增大至4.55 MPa;拱脚处环向应力较未冻胀环向应力增大约38.8%;拱底环向应力由0.5 MPa 增大至2.5 MPa;在衬砌外边界上,数值模拟所求环向应力与解析法所求环向应力吻合得较好,只是在α=100◦应力峰值附近有一定的偏差.

图5 衬砌外边界环向应力Fig.5 Circumferential stress along outer lining periphery

图6 为衬砌与冻胀圈接触面上的法向位移,从图中可以看到:数值法向位移解与解析法向位移解吻合得较好;衬砌受到冻胀力,产生了明显的冻胀变形.衬砌拱顶下沉较未受冻胀力时增大约37.4%.拱底隆起明显,与未受冻胀力时相比增大约3 mm.但是,总的来说计算所得未冻胀和冻胀法向位移均较实测值偏小[28],这可能与未考虑开挖卸荷因子和塑性变形有关.