Influences of ambient temperature on temperature field and mechanical behaviors of underground pipe galleries
-
摘要: 季冻区季节性温度变化对地下管廊的工程应用性能产生显著影响,为了研究季冻区环境温度对地下管廊温度场和力学特性的影响,采用ABAQUS有限元数值模拟技术分别研究了季节性环境温度作用下地下管廊温度场和土压力的变化规律。数值分析结果表明:管廊竖向截面最大冻深明显滞后于最低环境温度出现时刻;在升温与降温过程中,管廊上部及附近区域温度场先后出现“冷核区”和“热核区”等特殊现象,随着升温与降温的持续作用而逐渐消失;在降温与升温过程中,管廊顶板土压力先后出现增加与减小趋势;当季节性温度变化作用结束后,管廊顶板土压力显著大于初始状态。Abstract: Seasonal temperature changes in the seasonal frozen areas have a significant impact on the engineering application performance of the underground pipe galleries. To study the effects of ambient temperature in the seasonal frozen areas on the temperature field and mechanical behaviors of the underground pipe galleries, the change laws of temperature field and earth pressure of the underground pipeline galleries under seasonal ambient temperature are studied by conducting the ABAQUS finite element numerical simulation. The numerical analysis results indicate that the occurent moment of the maximum freezing depth of the vertical section inside the pipe gallery is lagged behind that of the lowest ambient temperature. During the process of heating and cooling, the temperature fields inside the upper and nearby areas of the pipe gallery successively exhibit the special phenomena such as "cold core zone" and "hot core zone". These phenomena gradually disappear with the continuous effects of heating and cooling. During the cooling and heating processes, the earth pressure on the roof of the pipe gallery shows an increasing and decreasing tendency. After the influences of seasonal temperature, the earth pressure on the roof is significantly greater than that at the initial state.
-
Keywords:
- ambient temperature /
- temperature field /
- mechanical behavior /
- pipe gallery /
- numerical simulation
-
0. 引言
地下洞室对弹性波的散射问题在地震工程、地下结构抗爆、地球物理勘探等领域有着广泛的应用。近年来,许多研究人员基于波函数展开法[1-2]对单相介质中衬砌结构的动力响应进行了探讨。然而,土和岩石具有多孔性和多相性,单相弹性连续介质无法准确描述实际情况。因此,波在两相介质中的传播和动力响应问题受到广泛关注。应用Biot饱和多孔介质波动理论,周香莲等[3]求解了全空间圆柱形衬砌洞室对平面波散射的解析解。Jiang等[4]引入大圆弧假设,解决了平面P波入射下饱和多孔弹性半空间中圆形衬砌洞室周围的散射问题。为了更符合实际工程,张海等[5]使用辅助函数配合波函数展开法和傅立叶级数展开方法,进一步研究了SH波入射下含直边界半圆形衬砌洞室的动力响应。Ding等[6]利用波函数展开法,给出了平面P波入射下饱和多孔弹性介质中复合衬砌动应力集中的解析解。Hasheminejad等[7]解析分析P波和SV波作用下饱和多孔弹性介质中变壁厚内衬的圆柱形洞室的动应力集中问题。Liu等[8]基于平面复变理论得到了多孔弹性半空间中衬砌洞室对P1和SV波散射的解析解。徐长节等[9],Xu等[10]基于非局部弹性理论分析了孔径效应和孔隙动力学效应对饱和土中圆柱形衬砌动力响应的影响。Ding等[11]基于非局部Biot理论得到了平面P波入射下饱和多孔介质中浅埋复合衬砌洞室动力响应的解析解。
综上,饱和多孔介质中衬砌洞室的动力响应已经得到了广泛而充分的研究。可实际上,土是一种三相介质,即使地下水位以下的土,也可能存在不完全饱和的情况。Li等[12]总结了众多非饱和孔隙弹性模型以及非饱和多孔介质波动特性的研究,表明由于气相的存在,非饱和介质的动态响应在许多方面与传统饱和多孔介质的动态反应显著不同。然而,目前有关非饱和多孔介质中衬砌洞室对地震波散射问题的研究还较为匮乏。Li等[13]研究了非饱和多孔介质中圆形洞室在内部荷载作用下的动力响应。Tan等[14]研究了无限非饱和多孔弹性介质中不完全接触的衬砌洞室对平面P波的散射问题。以考虑混合物理论的非饱和多孔介质动力理论为基础,采用波函数展开法和分离变量法,分别推导了平面P波和SV波入射时衬砌洞室周围位移应力的解析解。分析了不同入射频率下饱和度对衬砌洞室周围动应力集中系数的影响。
1. 分析模型
1.1 数学模型和基本方程
为了简化研究模型,把非饱和多孔介质中圆柱形衬砌洞室,视为无限非饱和多孔介质全空间中一个无限长度的衬砌洞室(图 1)。
衬砌洞室的内径为a,外径为b。本研究基于Wei等[15]提出的非饱和多孔介质动力学理论。控制方程表示为
nS0ρS0¨uS=(MSS+G)∇∇⋅uS+G∇⋅∇uS+MSW∇∇⋅uW+MSN∇∇⋅uN+ˆμW(˙uW−˙uS)+ˆμN(˙uN−˙uS),nW0ρW0¨uW=MSW∇∇⋅uS+MWW∇∇⋅uW+MWN∇∇⋅uN−ˆμW(˙uW−˙uS),nN0ρN0u ¨N=MSN∇∇⋅uS+MWN∇∇⋅uW+MNN∇∇⋅uN−ˆμN(˙uN−˙uS)。} (1) 式中:上角标S,W,N分别为非饱和介质中的固体骨架部分、液相部分和气相部分;G为非饱和多孔介质的剪切模量;nα0,ρα0,uα分别为各相的初始体积分数、初始密度和初始位移;MSS,MWW,MNN,MSW,MSN,MWN为与体积模量、剪切模量及模量间耦合相关的弹性系数;ˆμf为渗透率相关的系数,f=W,N。Wei等[15]给出了相关参数的表达式。采用Van等[16]提出的土水特征曲线V-G模型。
非饱和多孔介质中的本构关系表示为
σS=(MSS∇⋅uS+MSW∇⋅uW+MSN∇⋅uN)I+2G[∇uS+(∇uS)T],σw=(MSW∇⋅uS+MWW∇⋅uW+MWN∇⋅uN)I ,σN=(MSN∇⋅uS+MWN∇⋅uW+MNN∇⋅uN)I 。} (2) 式中:σα为各相的应力张量,α=S,W,N;I为单位张量矩阵。
衬砌介质被视为连续、各向同性的单相弹性介质,衬砌介质位移为uL,衬砌介质的拉梅常数为λL,μL,衬砌密度为ρL。
1.2 动力学方程的求解
根据亥姆霍兹(Helmholtz)矢量分解原理,非饱和多孔介质方程中的位移可以表示为两部分:
uα = ∇φα + ∇×Ψα(α=S,W,N)。 (3) 式中:φα,Ψα分别为非饱和多孔介质中α相的标量势函数和矢量势函数。
本文考虑简谐波入射,势函数可表示为φ= φ(x,y,z)e−iωt。将式(3)代入方程(1)可得
(∇2−k2j)φSj=0,φWj=δWjφSj,φNj=δNjφSj(j=1,2,3),} (4a) (∇2−k24)ΨS=0,ΨW=ζWΨS,ΨN=ζNΨS。} (4b) 从式(4a),(4b)中可以看出非饱和多孔介质中有3种压缩波1种剪切波,ki为3种压缩波波数, i=1,2,3,k4为剪切波的波数,δWj,δNj分别为液相、气相压缩波参与系数,j=1,2,3; ζW,ζN分别为液相、气相剪切波参与系数。详细表达式参见文献[13]。
将式(4a),(4b)中固相波函数转换到圆柱坐标系下并使用分离变量法求解,可得固体骨架的标量势函数和矢量势函数的通解为
φSj(r,θ)=∞∑n=0An,jKn(kjr)cosnθ (j=1,2,3), (5) ΨS(r,θ)=∞∑n=0BnKn(k4r)sinnθ。 (6) 式中:An,j,Bn为待定系数;Kn为第二类修正Bessel函数,下标n代表阶数。
同理可以得到衬砌介质中标量势函数和矢量势函数的通解:
φL(r,θ)=∞∑n=0(A(1)nIn(kL1r)+A(2)nKn(kL1r))einθ, (7) ΨL(r,θ)=∞∑n=0(B(1)nIn(kL2r)+B(2)nKn(kL2r))einθ。 (8) 式中:φL,ΨL分别表示衬砌介质中的标量势函数和矢量势函数;A(1)n,A(2)n,B(1)n和B(2)n为待定系数,分别表示P和SV波的振幅;In为第一类虚宗量Bessel函数,下标n代表阶数。
2. 波场分析与求解
2.1 波场分析
假设具有圆频率的平面谐波P1或SV波沿x轴正方向传播,用贝塞尔-傅立叶级数表示为
入射P1波:φi=φ0∞∑n=0εninJn(k0r)cos nθe−iωt, (9) 入射SV波:Ψi=Ψ0∞∑n=0εninJn(k0r)cos nθe−iωt。 (10) 式中:k0=ω/cp(入射P1波)或者k0=ω/cs(入射SV波)表示入射波的波数,φ0(Ψ0)表示入射波振幅,Jn(.)表示第一类n阶贝塞尔函数,若n=0则εn=1,若n⩾1则εn=2且i=√−1,t表示时间,e−iωt为时间因子将从后续表达式中省略。
入射平面波进入非饱和多孔介质,由于衬砌洞室的存在,会在非饱和多孔介质中产生散射波,在衬砌介质中产生透射波和内表面的反射波。以P1波入射为例,非饱和多孔介质中的波场可表示为
φS=φ0∞∑n=0εninJn(k0r)cos nθ+∞∑n=0An,jKn(kjr)cosnθ,φW=δW1φ0∞∑n=0εninJn(k0r)cos nθ+3∑j=1δWj∞∑n=0An,jKn(kjr)cosnθ,φN=δN1φ0∞∑n=0εninJn(k0r)cos nθ+3∑j=1δNj∞∑n=0An,jKn(kjr)cosnθ ,} (11a) ΨS=∞∑n=0BnKn(k4r)sinnθ,ΨW=ζWΨS,ΨN=ζNΨS。} (11b) 衬砌介质中的总波场可表示为
φL=∞∑n=0(A(1)nIn(kL1r)+A(2)nKn(kL1r))cosnθ,ΨL=∞∑n=0(B(1)nIn(kL2r)+B(2)nKn(kL2r))sinnθ。} (12) 2.2 边界条件
当r=a时,在衬砌介质中:
σLrr=0,σLrθ=0。} (13) 当r=b时,非饱和多孔介质与衬砌介质满足连续变形条件与界面不透水边界条件:
urL=urS=urW=urN,uθL=uθS,σrθL=σrθS,σrrL=σrrS+σW+σN。} (14) 式中:uαr,uαθ分别为非饱和多孔介质各相位移法向和切向分量,α=S,W,N;urL,uθL分别为衬砌的法向、切向位移分量;σrrS,σrθS分别为非饱和多孔介质固相法向、切向应力分量;σW,σN分别为非饱和多孔介质的孔隙水压力和孔隙气压力;σrrL,σrθL分别为衬砌法向和切向应力分量。
2.3 求解
将非饱和多孔介质和衬砌介质中总波场势函数的式(11),(12),代入应力和位移的式(2)中,根据边界条件,可得到波场中未知幅值系数的方程:
[G1G2G3G4G5G6G7G8H1H2H3H4H5H6H7H8T1T2T3T40000U1U2U3U40000V1V2V3V5−D1|r=b−D2|r=b−D3|r=b−D4|r=bY1Y2Y3Y5−F1|r=b−F2|r=b−F3|r=b−F4|r=b0000D1|r=aD2|r=aD3|r=aD4|r=a0000F1|r=aF2|r=aF3|r=aF4|r=a][An,1An,2An,3BnA(1)nA(2)nB(1)nB(2)n]=[G9H9T5U5−V4−Y400]. (15) 系数Gi,Hi,Di,Fi,Ti,Ui,Vi和Yi的表达式由上述公式计算得出,表达式略。通过求解式(15),得到未知波幅系数,以确定非饱和多孔介质和衬砌中的波场,从而获得其中位移和应力的解析解答。
3. 验证和算例
3.1 退化验证
将衬砌外内半径之比取1.000001,其结果将会退化为非饱和多孔介质中无衬砌洞室的情况。与文献[13]计算结果对比,取参数(n0=0.34,Sr=0.6),图 2给出了P波入射下衬砌洞室与空腔洞室动应力集中因子对比结果。从图 2中可以看出,当衬砌很薄时,衬砌洞室退化后的结果与洞室的计算结果一致。
3.2 算例分析
工程抗震研究的重点是衬砌的动应力集中因子(dynamic stress concentration factor,简称DSCF),本文中考虑介质中环向应力与入射波产生的环向应力的比值,表达式为σθθS/σ0,其中
σθθS=(MSS∇2φS+MSW∇2φW+MSN∇2φN)+2Gr(∂φS∂r+1r∂2φS∂θ2)−2G[∂∂r(1r∂ΨS∂θ)], (16) σ0=[2G+MSS+MSW+MSN+δW1(MWW+MSW+MWN)+δN1(MSN+MWN+MNN)]k21 (P波入射), (17a) σ0=Gk24(SV波入射)。 (17b) 从应力表达式中可以看到动剪切模量是一个关键因素,在以往的研究中很少将饱和度作为动剪切模量的考量,实际饱和度的变化对土中动剪切模量影响很大[17]。采用下式计算非饱和土动剪切模量[17]:
G=G0+2050αln(√(Se)−2−1+(Se)−1)tanϕ′。 (18) 式中:G0为土体饱和时动剪切模量;ϕ′为土体饱和时内摩擦角。
本文分别分析P波和SV波入射时,不同频率下饱和度对衬砌洞室周围动应力集中因子的影响。为避免孔洞尺寸对分析结果的影响采用无量纲频率η=a×k0,式中a为衬砌内半径,k0为入射波波数。衬砌介质的参数取λL=3.5×107 kPa,μL=3.5×107 kPa,ρL=2500 kg/m3。非饱和多孔介质的材料参数及土水特征曲线参数α,d,SrW ,SrN ,如表 1所示[13, 17]。
表 1 非饱和多孔介质材料参数Table 1. Material parameters of unsaturated porous medium孔隙率n0 ρS0/(kg·m-3) ρW0/(kg·m-3) ρN0/(kg·m-3) KS/kPa 0.36 2650 1000 1.1 3.6×107 KW/kPa KN/kPa νW/(Pa⋅s) νN/(Pa⋅s) G0/kPa 2×106 110 1.0×10-3 1.8×10-5 9.23×107 ˆK/kPa κ/(m⋅s−1) αB SrW SrN 2.0×105 2.5×10−12 1 0.05 1 υ d α ϕ′/(°) 0.3 2 2.0×10−5 20 取衬砌外内半径之比为ϑ=1.1,分别计算η为0.25,1.0,2.0情况下饱和度为Sr为0.2,0.6,0.8,1.0时衬砌外壁动应力集中因子分布情况。图 3,4分别给出了P波和SV波入射时衬砌洞室周边动力集中因子分布曲线。从图中可以看出,P波和SV波入射时衬砌周边的应力均对称分布在入射波两侧,饱和介质与非饱和介质中动应力分布有明显差异。P波入射时,饱和度变化对动应力分布形状影响较大,饱和时动应力分布形状与非饱和状态时有显著不同;SV波入射时,饱和度变化对动应力分布形状影响较小,饱和时动应力大于非饱和状态。P波入射下,当η=0.25时,Sr < 1.0动应力集中因子的最大值随饱和度的增加逐渐增大,在Sr为0.2,0.6时分布形式较一致呈“8”字形,最大值出现在80°,280°;在Sr=0.8时图像分布形式变为3个凸起,最大值出现在70°,180°和290°。Sr=1.0时动应力集中因子呈椭圆形沿洞中均匀分布,最大值角度出现在80°,280°。当η=1.0时,Sr < 1.0动应力集中因子最大值随饱和度的增加逐渐增大,动应力集中因子分布呈“豌豆状”,随饱和度增加逐渐饱满。Sr=1.0时动应力集中因子呈椭圆形分布在坐标系的中上部,最大值出现在0°,360°。当η=2.0时,动应力集中因子最大值随饱和度的增加逐渐增大,曲线集中分布在坐标系的上半部分,随饱和度增加最大值出现位置逐渐向0°坐标轴靠拢。P波入射时,最小值位置出现在入射面和入射背面,基本不受饱和度和入射频率的影响。
由于P波与SV波入射的振动形式不同,衬砌洞室周边动应力集中因子分布不同,在SV波入射时动应力集中因子呈蝶状分布。对比图 3,4可以发现不同相对频率下,饱和时动应力集中因子最大值要大于入射P波时的动应力集中因子。在图 4中可以发现,SV波入射时,饱和与非饱和情况存在明显差异,相比于Sr=1.0,Sr < 1.0时环向应力数值更小。在Sr < 1.0时饱和度对动应力集中因子影响较小,但随相对频率变化动应力集中因子变化明显,整体上随着入射频率的增加衬砌周边动应力逐渐增加,由于Sr=1.0与Sr < 1.0增幅不一致,在相对频率为2.0时饱和与非饱和动应力集中因子最接近。Sr < 1.0时最大值出现角度不变为40°,320°;Sr=1时最大值出现角度变化为50°,310°;30°,330°;20°,340°。SV波入射时最小值出现在轴线上。
4. 结论
推导了平面P波和SV波入射下非饱和多孔介质中内嵌圆柱形衬砌洞室的动力响应解析解。考虑了饱和度对剪切模量的影响,求解并分析了饱和度和频率对圆柱形衬砌洞室动应力集中因子的影响,得到3点结论。
(1)饱和度对动应力集中因子的影响不容忽视,无论P波还是SV波入射,饱和时动应力分布与非饱和时有显著不同。且P波入射时,饱和度对动应力集中因子的分布形式影响较大,随着饱和度增加,动应力因子最大值出现角度逐渐向0°偏移,动应力集中因子最大值逐渐增加;SV波入射下,饱和度小于1时,饱和度对动应力集中因子影响较小。
(2)由于P波与SV波的振动形式不同,衬砌洞室周边动应力集中因子分布不同。P波入射时,动应力集中因子分布更多呈椭圆形或“豌豆”状;SV波入射时,动应力集中因子分布呈蝶状。
(3)入射频率主要影响动应力集中因子的大小。随着入射频率的增加,P波和SV波入射下洞室周边动应力集中因子最大值均逐渐增加。
-
表 1 混凝土面层与管廊的力学与热物理参数
Table 1 Mechanical properties and thermal physical parameters of concrete surface and pipe gallery
参数 值 密度ρ/(g·cm-3) 2.4 弹性模量E/MPa 20000 泊松比 0.25 热传导系数λ/(W/(m·℃)) 1.5 比热容C/(J/(kg·℃)) 970 表 2 土体热膨胀系数
Table 2 Thermal expansion coefficients of soil
T/℃ -20 -10 -3 0 > 0 α/(℃-1) -6×10-4 -7.2×10-4 -9.6×10-4 -1.2×10-3 0 注:负号表示土体在低温作用下体积膨胀,0表示土体在正温作用下体积保持不变。 表 3 土体力学与热物理参数
Table 3 Mechanical and thermal physical parameters of soil
主要参数 值 冻土 融土 黏聚力c/kPa 12 12 内摩擦角ϕ/(°) 22 22 密度ρ/(g·cm-3) 1.89 1.89 弹性模量E/MPa 75 35 泊松比 0.25 0.3 热传导系数λ/(W/(m·℃)) 1.82 1.63 比热容C/(J/(kg·℃)) 2720 3140 -
[1] LAI J X, QIU J L, FAN H B, et al. Freeze-proof method and test verification of a cold region tunnel employing electric heat tracing[J]. Tunnelling and Underground Space Technology, 2016, 60: 56-65. doi: 10.1016/j.tust.2016.08.002
[2] TAN X J, CHEN W Z, YANG D S, et al. Study on the influence of airflow on the temperature of the surrounding rock in a cold region tunnel and its application to insulation layer design[J]. Applied Thermal Engineering, 2014, 67(1/2): 320-334.
[3] LAI Y M, WU Z W, ZHU Y L, et al. Nonlinear analysis for the coupled problem of temperature, seepage and stress fields in cold-region tunnels[J]. Tunnelling and Underground Space Technology, 1998, 13(4): 435-440. doi: 10.1016/S0886-7798(98)00086-8
[4] LIU L L, LI Z, LIU X Y, et al. Frost front research of a cold-region tunnel considering ventilation based on a physical model test[J]. Tunnelling and Underground Space Technology, 2018, 77: 261-279. doi: 10.1016/j.tust.2018.04.011
[5] XIAO C Z, CUI F L, DING L Q, et al. Temperature distributions in geogrid-reinforced soil retaining walls subjected to seasonal freeze-thaw cycles[J]. International Journal of Geomechanics, 2022, 22(12): 04022234. doi: 10.1061/(ASCE)GM.1943-5622.0002595
[6] YARIVAND A, BEHNIA C, BAKHTIYARI S, et al. Performance of geosynthetic reinforced soil bridge abutments with modular block facing under fire scenarios[J]. Computers and Geotechnics, 2017, 85: 28-40. doi: 10.1016/j.compgeo.2016.12.004
[7] 高煜, 刘勇, 石苏意, 等. 泡沫轻质土路堤稳定性分析及材料参数控制标准[J]. 长沙理工大学学报(自然科学版), 2017, 14(1): 24-30, 83. doi: 10.3969/j.issn.1672-9331.2017.01.004 GAO Yu, LIU Yong, SHI Suyi, et al. Foam light weight soil embankment stability analysis and control standard of material parameters[J]. Journal of Changsha University of Science and Technology(Natural Science), 2017, 14(1): 24-30, 83. (in Chinese) doi: 10.3969/j.issn.1672-9331.2017.01.004
[8] 中华人民共和国建设部. 砌体结构设计规范: GB 50003—2001[S]. 北京: 中国建筑工业出版社, 2002. Ministry of Construction of the People's Republic of China. Code for Design of Masonry Structures: GB 50003—2001[S]. Beijing: China Architecture & Building Press, 2002. (in Chinese)
[9] 付晓丹. 季冻区悬臂式支挡路基温度场与应力场数值分析[D]. 兰州: 兰州交通大学, 2019. FU Xiaodan. Numerical Analysis of Temperature and Stress Fields in Subgrade with Cantilever Retaining wall in Seasonally Frozen Regions[D]. Lanzhou: Lanzhou Jiaotong University, 2019. (in Chinese)
[10] LAI Y, WANG Q, NIU F, et al. Three-dimensional nonlinear analysis for temperature characteristic of ventilated embankment in permafrost regions[J]. Cold Regions Science and Technology, 2004, 38(2/3): 165-184.