寒区引水隧洞围岩冻胀特性分析——以新疆布伦口水电站引水隧洞工程为例

彭小丽，王海娟，杨丰春，姜海波，李唱唱

(石河子大学水利建筑工程学院，新疆 石河子 832003)

0 引言

地球上冻土区面积约占陆地面积的50%，我国多年冻土和季节冻土面积分别占国土面积的21.5%和53.5%[1]。随着“一带一路”倡议深入实施，越来越多的隧洞需要建立在我国西部寒区，因此冻岩问题已经成为近年来的研究热点。

寒区隧洞普遍受低温、渗流共同作用而产生冻胀力。低温冻融环境下，围岩因温度梯度进行热传导而使得内部围岩温度场发生改变，围岩裂隙中的液态水状态也发生改变。当温度达到冰点时，孔隙水冻结，低温冻结作用会产生9%的体积膨胀。杨更社等[2]对软岩材料进行了在一定温度梯度下的水热试验，结果表明温度梯度是岩体水分迁移的主要动力，岩石孔隙率与冻结时间成正相关。谭贤君等[3]建立含相变的水热耦合模型，研究了隧道温度场分布规律和冻融圈大小。夏才初等[4]基于稳态假定，积分推导了考虑围岩、衬砌未冻水含量的最大冻结深度解析计算式，并通过温度场反演获取解析式中的参数从而求得最大冻结深度半解析解。

在冻土和冻岩的冻胀特性研究方面，国内外许多学者采用室内试验、数值模拟和现场监测等方法相结合，取得了丰硕的成果。J.M.Konrad[5]提出了冻土分凝势理论，并根据分凝势理论研究了冻结方式对冻结特性的影响。赖远明等[6]利用黏弹性原理推导出冻胀力求解相关公式，并采用数值逆变换方法求出寒区隧洞冻胀力，通过算例计算表明是否考虑冻胀因素对应力值大小存在着较大的差异。张德华等[7]、Gao等[8]、黄继辉等[9]基于弹塑性理论，给出了不同应力准则下围岩冻胀力的弹性力学解析解。康永水等[10]运用应变片法对低温环境下饱和岩样和干燥岩样进行测试，得到岩样的低温应变特征，研究其冻胀融缩效应，并给出岩石冻胀变形规律。渠孟飞等[11]、严健等[12]为研究冻胀压力和冻融圈厚度变化，分别进行了室内隧道模型试验和现场原位测试试验。耿珂[13]、吴剑等[14]依靠数值模拟方法对隧道围岩冻胀力和冻胀特性进行研究。

寒区隧洞围岩因水-冰相变的影响，在环境温度、孔隙渗流以及应力三者之间相互影响、相互作用下表现出“冻胀”现象，因冻胀作用而产生的隧洞冻害严重时会影响隧洞的施工和运营。大多数学者在冻胀数值模拟中仅考虑了热力耦合的过程，而未将孔隙水相变过程考虑在内。鉴于此，本文以新疆某水电站引水隧洞为研究对象，使用大型有限元模拟软件，考虑水-冰相变现象，建立低温冻结环境下的围岩瞬态温度传导模型和温度-渗流-应力三场耦合模型，计算冻结深度和冻胀力特性，并将冻胀发生前后的冻胀特性进行对比，以此分析冻胀作用对该引水隧洞的影响，以期为引水隧洞的安全施工和运营提供一定指导。

1 冻胀力弹塑性力学解析解

关于围岩冻胀模型，存在多种研究理论[15]。本文采用围岩整体冻胀模型，认为受温度影响围岩沿径向范围内形成冻融圈，当温度降低时，围岩孔隙中的水结冰使得冻融圈整体膨胀，因此产生了冻胀力。

寒区隧洞冻胀力弹塑性力学解析解求解模型如图1所示。其中，ph为围岩冻胀后冻结层内壁作用在衬砌上的围岩压力;pf为冻胀后冻结层外壁作用在未冻围岩上的围岩压力;hf为冻结深度;r0、rⅠ、rⅡ分别为衬砌内径、衬砌外径(或冻结层内径)、冻结层外径;rⅢ为任意一处未冻结围岩距圆心的距离。为简化计算并使解析解的值更接近实际值，需对模型作出如下假设： 1)假定衬砌、已冻结岩、未冻结岩均为弹性体； 2)衬砌与已冻结围岩为绑定接触，冻融圈与原始围岩为弹性接触； 3)各材料符合各向同性假设； 4)不考虑风化破碎层以及孔隙水渗流的耦合作用。

图1 弹性力学求解模型

区域Ⅰ为圆形隧洞衬砌，将衬砌看作受外压ph作用的圆筒结构，根据弹性力学[16]，该区域应力场拉密解可表示为

(1)

径向位移

(2)

式中EⅠ、μⅠ分别为衬砌的弹性模量和泊松比。

当r=rⅠ时，衬砌外壁位移

(3)

区域Ⅱ为已冻结围岩，已冻结区受ph、pf共同作用，该区域内应力拉密解为

(4)

其径向位移

(5)

式中：rⅡ=rⅠ+hf；hf为冻结深度，m，计算时hf的取值通过对拱腰位置布置温度测量元件进行围岩内部径向温度场监测，确定为2.35 m；EⅡ、μⅡ分别为已冻结岩的弹性模量和泊松比。

当r=rⅠ时，冻结圈内壁发生的位移

(6)

当r=rⅡ时，冻结圈外壁发生的位移

(7)

区域Ⅲ为未冻结围岩，未冻围岩外半径趋于无穷大时，令无穷远处冻胀压力p=0，该区域内的拉密解为

(8)

径向位移

(9)

式中EⅢ、μⅢ分别为未冻结围岩的弹性模量和泊松比。

当r=rⅡ时，未冻结围岩内壁发生的位移

(10)

根据位移连续条件，衬砌与冻融圈内壁(r=rⅠ)、冻融圈外壁与未冻结围岩内壁(r=rⅡ)需满足的位移条件为[17]：

(11)

式中Δ1、Δ2分别为冻胀后已冻结围岩内壁和外壁的变形量，表达式见式(2)。

(12)

式中ΔV为围岩冻胀前后的体积差。

由几何关系，径向冻胀变形量可表示为

(13)

式中ηr为考虑冻结过程中岩石自身约束作用的饱和冻胀率[18]，表达式见式(14)。

ηr=2.17%n。

(14)

式中n为围岩孔隙率。

联立式(1)—(14)，求得冻结围岩的冻胀力弹性力学解析解为：

(15)

(16)

本文考虑到围岩和围岩孔隙水在冻胀受冻胀过程中岩体骨架的约束作用，对径向冻胀变形量式(13)中的冻胀率进行修正，并代入其他工程数据对修正后的冻胀力解析式进行对比以验证公式的正确性，如表1所示。

表1 其他冻胀力解析计算值与本文解析计算值对比

通过将如表1所示的4个文献中的参数代入到本文公式中进行计算，并与其他文献中的计算值及实测值进行对比，得出利用本文所推导的公式计算出的冻胀力理论值接近其实测值。

利用式(15)—(16)，对新疆布伦口水电站引水隧洞进口50 m处断面进行计算，得到冻结圈围岩作用在衬砌上的冻胀力ph=10.856 MPa，作用在未冻结围岩上的冻胀力为pf=4.002 MPa。将所求ph和pf值代回到式(5)、式(6)中，即可得到冻结圈围岩受冻胀作用应力场的σr和σθ分别为：

(17)

2 基于水-冰相变的引水隧洞围岩温度场分析

2.1 工程地质概况

新疆布伦口水电站引水隧洞工程全长20.14 km，地质情况较为复杂。根据地质测绘，坝址区发育规模较大的断裂主要为F1、F2、F5。引水隧洞位于高寒半干旱气候区，季节性温差较大，雨季和旱季明显，年降水量不大，洞址区地下水多以基岩裂隙水为主。隧洞洞址处多年平均气温为0.7 ℃，年最低温度为-34.3 ℃，年最高温度为 35.9 ℃，最冷月平均气温为-13.5 ℃；多年平均降水量为127.5 mm，年平均蒸发量为2 297.9 mm；最大积雪厚度约160 mm，最大冻结深度为2.35 m。现场监测结果表明，沿隧洞轴线方向共计出现120多处结冰点，且其冻结程度不同。隧洞平导进口段部分断面衬砌漏水产生冻胀，洞内冻结现象严重，并且冻结的时间较长，影响隧洞的施工和运营，不利于隧道的结构安全。寒区隧洞是否发生冻胀破坏的首要因素是隧洞所处的温度环境，研究隧洞进口段的温度场和围岩冻胀力，是采取防冻害措施的重要依据。

2.2 相变模拟方法

寒区隧洞由于环境温度与围岩自身温度之间存在一定的温差，从而随时间产生瞬态温度传导，加之隧洞围岩中裂隙等结构面中水的存在会伴随温度的降低而发生相变。水冰相变影响热传导参数，设定水的相变域为(-0.5 ℃，0.5 ℃)，使用显热容法，将水-冰相变潜热的影响等效为围岩比热容的变化，可构造等价比热容和热传导系数[19]如下：

(18)

(19)

式(18)—(19)中：Cf、Cu分别为围岩中水融化和冻结时围岩比热容；λf、λu分别为融化和冻结时热传导系数；L为相变潜热；t为围岩温度。

2.3 有限元计算模型

引水隧洞洞径D=5 m，越靠近进水口其围岩破碎程度越高。对该引水隧洞现场温度监测分析可知，越往隧洞内部越不易发生冻胀[20]，故选取隧洞进口段洞深50 m处的截面作为计算断面，对围岩温度场200 d的变化情况进行分析。

本文采用有限元模拟软件先对温度场传热进行模拟分析。在建立数值模型时根据圆形截面隧洞的对称性，选用隧洞断面的一半作为计算区域，围岩计算范围为22.5 m×40 m，断面总共划分748个单元，共计803个节点，有限元计算模型网格划分及尺寸如图2所示，图中a、b、c3点分别为拱顶、拱腰、拱底数据提取处。温度场传热过程为瞬态传热，计算总时长为200 d，每4 h设置1个增量步，共计800增量步。

(a) 模型网格 (b) 模型尺寸(单位： m)

根据工程地质资料及现场数据，考虑低温相变，围岩初始温度场为3.5 ℃，AE边界及隧洞洞腔BC边界的温度取10月首次出现负温至次年4月期间200 d内各月的平均温度，隧洞进口段月平均温度见表2，热流边界DF的热流密度为0.06 W/m2，其余热物理学参数见表3。

表2 隧洞进口段月平均气温

表3 围岩热物理参数

2.4 温度场结果分析

瞬态温度传导200 d后隧洞围岩温度发展规律如图3和图4所示。由图3可知，隧洞围岩温度场受洞腔内低温环境影响显著，围岩孔隙水受冻结发生相变，冻结锋面随着冻结时间的推移向围岩内部逐渐推移，表现为冻结区逐渐扩大，未冻区逐渐缩小。通过路径b—d(见图2)，从温度场分布云图中提取冻结深度结果，观察冻结深度随时间的变化规律(见图4)，自传热开始到传热时间为130 d时，径向冻结深度随冻结时间的增长而增加，最大冻结深度达2.04 m，与实际测得的冻结深度2.35 m存在一定误差，这是由于采用有限元数值模拟水冰相变瞬态传热过程时将围岩默认为各向同性介质，未考虑围岩内部实际存在的裂隙、节理结构，热传导系数会因结构面的存在而受到一定影响。随着传热过程进行至130 d左右，由于环境温度的升高，隧洞边界同冻结锋面温差Δt减小，加之已冻区和未冻区围岩导热系数不同，冻结面开始往外部移动，围岩已冻区逐渐减小，冻结深度开始减小。至180 d时，边界出现正温，围岩导热系数增大，温度变化加快，围岩温度全部变为正温区。

(a) 50 d (b) 100 d

(c) 150 d (d) 200 d

图4 冻结深度随时间的变化曲线

3 基于水-冰相变的围岩冻胀特性分析

岩土体产生冻胀需要满足以下3个条件： 1)冻胀敏感岩石； 2)冻结温度和时间； 3)初始水分或补给水分。因此，在计算冻胀力时温度传导和孔隙渗流也需要考虑在内。考虑到工程实际，围岩孔隙水在低温冻结条件下发生相变产生冻胀，冻结锋面在围岩内部移动，并且隧洞冻胀实际边界条件复杂，是一种强非线性问题。式(17)求得的冻胀作用下的应力场σr和σθ，可通过坐标转换求得σx和σy，进而解得主应力。但基于一定假设基础的理论解计算存在一定误差，且无法直接反映出围岩的应力分布状态。为反映围岩的应力状态并求得更接近于真实值的近似解，采用有限元软件对引水隧洞进行冻胀模拟，分析围岩冻胀特性。

3.1 计算模型及参数

建立考虑低温水冰相变的水热力三场耦合的冻胀特性计算模型(见图2)，对模型左右边界进行水平约束，下部边界进行竖向约束，上边界无约束，并在此基础上添加力学和渗流参数，进行多场耦合计算。但由于本文使用的有限元分析软件不能建立直接的温度-渗流-应力耦合关系，故采用温度场对渗流-应力耦合的影响来建立间接瞬态耦合模型进行冻胀计算。考虑低温水-冰相变，水的相变域为(-0.5 ℃，0.5 ℃)，相变潜热Lf=333.56 kJ/kg，热膨胀系数为8×10-6。假定岩体为弹塑性材料，服从Mohr-Coulumb屈服准则，根据相关规范和现场资料，计算所需的其他相关参数如表4所示。

表4 隧洞围岩力学计算参数

3.2 计算结果分析

围岩受冻胀200 d过程中的应力及位移变形如图5和图6所示，从数值分析软件中提取受冻胀前后隧洞拱顶、拱腰、拱底的应力及位移，对比变化如表3所示。观察图5可知，隧洞受冻胀后，拱腰周围围岩出现应力最大值，200 d内最大拱腰应力达0.635 MPa，拱顶和拱底位置应力较小。对比冻胀前后各部位应力值，拱顶应力由冻胀前的0.085 MPa增大到0.392 MPa，拱腰应力从0.604 MPa减小到0.527 MPa，拱底应力由0.115 MPa增加到0.399 MPa。冻胀后围岩拱顶、拱底应力值呈增大趋势，拱腰呈减小趋势，各部位变化情况见表5。

(a) 50 d (b) 100 d (c) 150 d (d) 200 d

(a) 冻胀前水平位移 (b) 200 d水平位移 (c) 冻胀前竖直位移 (d) 200 d竖直位移

表5 冻胀前后围岩应力及变形对比

观察图6(a)—(b)可知，受冻胀影响，隧洞拱顶和拱底处水平位移因受到水平方向的位移约束无明显变化；从有限元分析软件后处理中提取冻胀前后隧洞拱腰周围的水平位移，呈水平收缩趋势，变为-0.36 mm。由图6(c)—(d)可知，围岩拱顶竖直方向上位移发生沉降，位移由冻胀前的-9.59 mm变为冻胀后的-8.41 mm，竖向位移变化量为1.18 mm。

综上，受冻胀影响，围岩拱顶、拱底应力呈增大的趋势，拱腰应力逐渐减小，洞壁变形呈收缩趋势。随冻结时间的增加，冻胀引发的冻深变化、应力变化和位移变形会进一步使得围岩结构发生冻害，进而威胁围岩稳定性。

4 结论与讨论

利用弹性力学理论对围岩冻胀力解析解进行计算，基于有限元分析软件，考虑水-冰相变下的温度场、渗流场、应力场相互作用，建立有限元模型对冻胀特性进行数值模拟分析计算，得出如下结论：

1)引水隧洞围岩在长时间的低温作用下，孔隙水冻结发生相变，随时间的推移，围岩冻结区增大，冻结深度随之增加，低温冻结200 d内最大冻结深度为2.04 m。由于冻结深度变化速度受围岩温度梯度及环境温度影响，130 d后冻结深度逐渐减小直至全部变为正温区。

2)在相变引起的冻胀作用下，拱顶应力由冻胀前的0.085 MPa增大到0.392 MPa，拱腰应力从0.604 MPa减小到0.527 MPa，拱底应力由0.115 MPa增加到0.399 MPa，对比冻结前后应力状态，拱顶、拱腰、拱底的变化幅度分别为361.2%、-12.7%、247.0%。拱腰处水平位移变为-0.36 mm，表现为水平收缩，围岩拱顶竖直方向发生沉降，产生1.18 mm变化量。

寒区水工隧洞在低温水-冰相变影响下冻深逐渐增大，围岩的应力和位移受低温冻结作用影响发生明显改变，为保证围岩稳定性和隧洞的安全运营，应及时对隧洞施加合理的保温和支护措施。因此，在本文对围岩冻胀特性研究基础上，后续可通过数值模拟、现场监测等方法对保温层和衬砌进行进一步研究。