水热力耦合下公路隧道结构各向异性冻胀力分布特征研究

陆许峰， 黄解放， 张 锋， 冯德成

（1. 哈尔滨工业大学 交通科学与工程学院，黑龙江 哈尔滨 150090； 2. 中交一公局集团有限公司，北京 100024）

0 引言

随着我国东北及西部地区交通事业的发展，大量地处寒区的隧道得以修建。然而，寒区隧道的冻胀力给衬砌结构带来了极大挑战。明确隧道冻胀力的分布，将其更好地考虑至隧道结构设计中，是提高隧道抗冻耐久性的关键。

现有的隧道冻胀模型大体可分为含水风化层冻胀模型、衬砌背后局部积水冻胀模型，以及冻融圈层整体冻胀模型［1-2］。其中含水风化层冻胀模型主要关注拱肩至边墙10～20 cm 范围内风化层的冻胀［3］，一般将洞口段与洞身段的冻胀力分别取为0.9 MPa 与0.6 MPa［3-4］。然而，风化层本身是围岩冻胀导致的结果而非诱因［5］，且冻胀力实际为冻深范围内的风化层与扰动层共同产生［6-8］，该类模型依旧有一定的提升空间。局部积水冻胀模型则假定衬砌背后有存水空间，并将存水空间假定为一定的几何形状［9-10］，以积水空间内水膨胀诱发的力作为局部荷载。对于土质围岩、软岩或较破碎岩体等低模量工况，含水风化层冻胀模型只考虑了衬砌两侧风化层的冻胀，局部积水冻胀模型计算出的冻胀力较小，需要从冻融圈层整体冻胀的角度对衬砌-已冻围岩-未冻围岩系统进行分析［1，11-12］。冻融圈层整体冻胀模型大多在极坐标系下进行受力分析［13-14］，且衬砌与冻结圈大多被简化为圆形［11，15］。考虑非圆形衬砌几何，可使隧道冻胀力分布结果更为符合实际［16-17］。

隧道冻胀力的分布在衬砌外围一般是不均匀的。围岩单元在空间上的冻胀变形呈现出正交各向异性［18-19］。现有研究结果表明不均匀冻胀系数的增大会使冻胀力增大［20-21］，且考虑不均匀冻胀系数后的冻胀力计算结果也更为符合实际［22］。然而，即使考虑了衬砌几何形状与围岩冻胀的各向异性，从纯力学角度指定冻结圈来计算冻胀力，依旧无法完全解释隧道的冻胀力分布［2］。事实上，围岩冻结过程中的水-热-力耦合效应普遍存在。赖远明等［23］认为应该充分考虑渗流场与应力场对冻结深度的影响。Feng等［24］通过考虑围压与轴压的冻胀试验研究发现，应力边界对土体冻胀变形有着显著的影响。耿珂等［25］开展了内蒙古阿拉坦公路隧道运营期冻胀力的热-力耦合数值模拟，但未考虑冻胀变形的各向异性。Li 等［26］构建了隧道结构的水-热-力耦合模型，发现边墙处的积水导致该处的冻胀力更大。张泽等［27］发现隧道围岩水-热耦合过程中的相变潜热对温度场影响较大。王志杰等［28］发现围岩相变潜热对衬砌结构的主应力分布与大小均存在影响。

综上可知，冻融圈层整体冻胀模型相较其他两类模型而言具备更高的适用性。然而，隧道冻胀力的分布受到冻胀变形的各向异性、隧道几何，以及多场耦合效应的影响［18，23，29］，呈现出一定的复杂性。一方面，隧道围岩冻胀过程中的水-热-力耦合效应客观存在，在研究中不容忽视；另一方面，岩土体的冻胀同时发生在温度梯度方向与垂直于温度梯度方向，将体冻胀率分配至三维空间中更为符合实际，与公路等层状结构不同，这在隧道这类近似环状的结构中更需要得到重视。针对隧道冻胀力的分布，目前尚未能有在水-热-力耦合框架下引入围岩各向异性冻胀变形的研究。理论模型的推导，以及数值模型的实现，均具有较大的难度。为了更好地分析并阐明隧道冻胀力的分布，本文首先构建了考虑正交各向异性冻胀变形的围岩冻胀水-热-力耦合理论模型，并依托COMSOL Multiphysics 有限元平台进行二次开发实现了理论模型的数值求解，验证了模型的可行性；进而，以东天山隧道为例，计算并分析了隧道的温度场、水分场与冻胀力分布。该研究对阐明季节冻土区隧道工程的冻胀力分布特征，揭示冻胀力差异分布的机理有着重要的理论意义与实际意义。

1 围岩冻胀水-热-力耦合理论模型

1.1 模型假设

寒冷地区的隧道温度场随着大气温度发生波动，且伴随着水分迁移与相变的发生。寒区隧道开挖后一般先喷射混凝土形成初期支护，之后进行二次衬砌的修筑，此类施工条件下衬砌与围岩的结合较为紧密。本研究重点关注来自围岩的冻胀力，面向围岩风化破碎程度较高、含泥质充填、呈石夹土状或土夹石状时的工况，参考相关研究［21，27，30-32］，做如下假设：

①围岩内部的水分以液态水的形式迁移，且符合达西定律；

②忽略冻结或融化过程中围岩骨架的热胀冷缩，仅考虑水冰相变的体积变形；

③忽略初始地应力、衬砌与围岩自重的影响，衬砌与围岩之间紧密结合；

④围岩的冻胀变形为各向异性，其余模型参数为各向同性，且模型符合平面应变假设。

1.2 隧道围岩的水-热耦合控制方程

根据能量守恒与质量守恒基本原理，围岩的水-热耦合方程如式（1）所示［27］：

式中：Cv为考虑相变潜热的等效体积热容（kJ·m-3·℃-1）；λ为导热系数（W·m-1·℃-1）；T为温度（℃）；L为水冰相变潜热（kJ·m-3）；ρi、ρw分别为冰、水的密度（kg·m-3）；θw、θi分别为体积含水量（%）与体积含冰量（%）；k为围岩渗透系数（m·s-1）；φm为围岩基质吸力水头（m）。

1.3 考虑围岩正交各向异性冻胀变形的应力-应变控制方程

隧道围岩冻结过程中，在垂直于热流的方向形成冻结锋面；随着冻结温度的渗入，围岩在热流方向产生的冻胀变形远大于热流正交方向。因此，围岩冻结呈现正交各向异性冻胀变形的特征：即热流方向是一个主方向，正交于热流的方向是另一个主方向，具体如图1所示。

图1 考虑围岩正交各向异性冻胀变形的坐标系示意Fig. 1 Schematic diagram of the coordinate system considering orthotropic frost heave deformation

在ab坐标系下，考虑围岩的正交各向异性冻胀变形，将冻胀应变分解为沿着热流梯度方向与垂直于热流方向的分量，应力-应变关系如式（2）所示。同时，根据平面应变假设，根据关系式εaz=εbz=εz=0，结合式（2），得到的ab坐标系下的应力-应变关系，可提取出其中的冻胀应变分量，具体如式（3）所示。

将式（3）中的冻胀应变分量按式（4）旋转至全局坐标系，可得到全局坐标系下的冻胀应变分量如式（5）所示。最终，全局坐标系下考虑围岩正交各向异性冻胀变形分量的围岩本构方程如式（6）所示。

式中：{εij}与{σij}分别为ab坐标系下的弹性矩阵的总应变分量与总应力分量；{}与{εij}el分别为冻胀应变分量与弹性应变分量；E、μ分别为围岩弹性模量与泊松比，其中冻结与未冻结状态下弹性模量分别取Ef与Eu；ξ为围岩冻胀变形的正交各向异性系数；为围岩的体冻胀率；θ为热流方向与全局坐标系x轴的夹角；[D]为弹性矩阵。

1.4 模型耦合参数

围岩冻胀过程中的水热力耦合动态参数［27］如式（7）所示：

式中：λs、λw、λi分别为围岩、水、冰的导热系数（W·m-1·℃-1）；θs为围岩的体积含量，θs= 1 -θw-θi；Cs为围岩骨架导热系数 （kJ·kg-1·℃-1），冻结状态下取Csf，未冻状态下取Csu；ks为饱和渗透系数 （m·s-1）；S为相对饱和度；θn、θm分别为残余体积含水率与饱和体积含水率；m、n为VG 模型参数，m=1-1/n；a0为进气值（m-1）；B(T)为固液比函数，可按式（8）计算［30］；ρ为围岩整体密度（kg·m-3）。

式中：a为固液比函数参数，取值与围岩破碎程度与碎块体间的填充土质相关。按徐敩祖等［33］的试验数据，砂土取值0.61，粉土取0.47，黏土取0.56。当土质中含有部分破碎岩石时，围岩的结冰量会有所降低，a值可结合实际工况略作折减。

冻胀体积应变与含冰量密切相关［30，34］，具体如式（9）所示：

式中：α、β为体冻胀率预估函数的拟合参数。

2 模型验证

2.1 模型几何与参数

为验证模型的可靠性，针对青沙山隧道的ZK33+970 断面进行相应的水热力耦合计算，该断面为V级围岩［22］。在COMSOL Multiphysics 有限元平台中进行几何输入与网格划分。为提高计算精度，衬砌外围冻融圈内网格进行了适当的加密，隧道网格划分如图2所示。

图2 青沙山隧道ZK33+970断面几何Fig. 2 Model geometry of Qingshashan tunnel at ZK33+970 section

参考文献［27，30-32］，模型参数如表1所示。衬砌仅导热无相变，体积热容为2 425 kJ·m-3·℃-1，导热系数为1.85 W·m-1·℃-1，弹性模量为28.5 GPa，泊松比为0.20。

表1 青沙山隧道模型围岩计算参数Table 1 Calculation parameters of surrounding rock of Qingshashan tunnel model

2.2 初始值与边界条件

隧道尺寸相较地层尺寸而言较小，结合青沙山隧道工况，选定初始地层体积含水率为28%，初始地层温度为2.9 ℃。同时本研究重点关注围岩对衬砌结构的冻胀力，忽略开挖、重力场与初始地应力等因素的影响，因此初始应力场与初始位移场均设置为0。

隧道冻融影响区域有限，将模型四周在导热与渗水方面均设置为零通量边界，同时法向位移为完全约束。大气温度与隧道衬砌内部温度边界一般可采用正弦函数［27，30］，将大气温度作为隧道内侧温度边界可反映近洞口段的实际工况，更好地反映季节冻土区隧道冻胀病害。衬砌内边界的温度结合青沙山隧道所在地实际气温与冻深设置为三角函数逐年波动的形式，温度均值为5.9 ℃，温度振幅为15 ℃，且衬砌内侧设置为自由位移边界。围岩与衬砌交界面为渗流物理场计算区域的内边界，无衬砌内部补水，设置为渗流零通量。

2.3 计算结果对比

隧道冻胀的水-热-力耦合模型主要依托COMSOL Multiphysics 有限元平台求解。采用2 套PDE模块分别表述渗流方程与热传导方程，可实现模型中的围岩水-热耦合计算。需要注意的是，衬砌区域仅进行热传导计算，围岩区域则同时进行渗流与热传导计算。水-热耦合计算可获得不同时刻的围岩含冰量分布，结合式（5）和式（9），可获得围岩冻胀应变分量。将冻胀应变分量输入固体力学模块的热膨胀接口，可实现隧道冻胀的水-热-力耦合数值计算。模型求解采用瞬态求解器，时间步进方法为向后差分，计算步长为10 d，共7 300 d。

图3为最大冻深时青沙山隧道冻胀力计算结果及其与现场数据［22］的对比。可以看出：（1）青沙山隧道整体的冻胀力分布为拱脚＞顶拱＞仰拱，最大冻胀力位于顶拱与拱脚过渡段，约为0.15 MPa，最小冻胀力位于仰拱中部；（2）冻胀力计算结果与现场监测整体吻合较好，整体上可反映青沙山隧道的冻胀力分布，也验证了水-热-力耦合隧道冻胀力模型的正确性。

图3 青沙山隧道冻胀力计算结果Fig. 3 Calculation results of frost heave force in Qingshashan tunnel

3 隧道结构冻胀力分布及特征

3.1 工程概况

东天山隧道位于新疆东部，是G575 线巴里坤-哈密公路的重要工程段。其中的ZK9+000～ZK9+320、K9+000～K9+310 段为Ⅴ级围岩。洞身浅埋段为F1 全新世活动断裂带及影响破碎带。围岩主要为晚更新统洪积粉质黏土，下伏基岩为凝灰质砂岩，受断层影响岩体极破碎，围岩稳定性差。隧址所在的哈密地区位于中纬度亚欧大陆腹地，为温带大陆性干旱气候。年平均温度为9.9 ℃，年平均降雨量为34.9 mm，年平均蒸发量为2 799.8 mm，年平均日照时间为3 591.2 h，无霜期198 d。

3.2 东天山隧道水热力耦合模型构建

3.2.1 隧道几何

东天山隧道主洞内轮廓采用受力条件好的单心圆形式。顶拱段的圆弧半径为555 cm，仰拱半径为1 500 cm，拱脚处采用半径为100 cm 的小曲率圆弧连接顶拱与仰拱。此外，初期支护为26 cm 厚C25喷射混凝土、二次衬砌为50 cm厚C30钢筋混凝土。东天山隧道所经区域为季节性冻土区域，除受到自上而下的大气对地层的冻结作用以外，隧道近洞口段一定范围内围岩热量散失最为明显。从距离隧道主洞洞口不同纵深断面的温度监测结果来看，在距洞口3 600 m 范围内的断面会受到大气负温的影响。以距离洞口191 m纵深处的ZK9+000截面为例，选取隧道中心上下各60 m、左右各70 m 的范围作为计算区域。考虑到隧道横断面沿纵向延伸较长，将其简化为二维平面问题。图4 为东天山隧道主洞衬砌内轮廓几何，图5 为数值模型在隧道区域的网格剖分。

图4 东天山隧道主洞衬砌内轮廓尺寸图Fig. 4 The inner outline of Dongtianshan tunnel

图5 隧道区域网格剖分Fig. 5 Meshing of the tunnel model

图6 新疆哈密地区近20年气温拟合结果Fig. 6 Fitting results of temperatures in Hami area， Xinjiang in the past 20 years

3.2.2 初始值与边界条件

本文重点关注冻融圈层范围内围岩冻胀的影响。东天山隧道衬砌的最大半径为15 m，且季节冻土区隧道的冻融圈层深度也一般不超过2 m［22］。冻融区域相较地层尺寸较小。为简化计算将地层初始值做简单均匀化处理。结合现场调查结果，围岩的初始体积含水率取30%，初始温度取5 ℃。计算截面距离洞口较近，将大气温度作为隧道近洞口段内部温度边界，可更好地揭示东天山隧道的冻胀分布特征。衬砌内部的温度边界采用新疆哈密地区近20 年日均温拟合结果。模型的四周外边界为导热、渗水零通量，且法向位移完全约束。隧道内部的力学边界设置为自由。

式中：T为新疆哈密地区日均温（℃）；t为从每年1月1日开始算起的天数。

3.2.3 参数取值

结合现场资料与现有研究成果［27，30-32］，围岩的模型计算参数如表2 所示。围岩服从摩尔-库伦强度屈服准则。初期支护的导热系数与体积热容分别为1.80 W·m-1·℃-1与2 300 kJ·m-3·℃-1，弹性模量为20 GPa，泊松比为0.20；二次衬砌的导热系数与体积热容分别为1.85 W·m-1·℃-1与2 425 kJ·m-3·℃-1，弹性模量为31 GPa，泊松比为0.20。

表2 东天山隧道模型围岩计算参数Table 2 Calculation parameters of surrounding rock of Dongtianshan tunnel model

3.2.4 计算结果

图7 为东天山隧道运营20 年（7 300 d）内不同深度处的温度分布。可以看出：（1）在大气温度的影响下，衬砌及围岩的温度均发生着周期性的波动，且波动幅度随着深度的增加而逐渐减小。（2）在距离二衬内表面1 m 深度的围岩处，基本处于正温下的波动，但受负温轻度影响。距离二衬内表面3 m 处的围岩均处于正温下的波动。（3）距离二衬内表面1 m 深度范围内的围岩处于冻害影响区，若遇冻胀较为敏感、水分补给丰富的工况，有必要采取防排水、增设保温板等措施以提高隧道整体的抗冻耐久性。

图7 拱顶不同深度处的温度变化Fig. 7 Temperature changes at different depths of lining top

图8 为隧道不同位置处的体积含水率时程变化。可以看出：（1）随着衬砌内部大气温度的周期性波动变化，隧道不同位置围岩的含水率均发生着周期性的波动。（2）单个冻结周期内，围岩的体积含水率变化可分为冻结、融化、滞水、滞水消散4 个阶段。（3）随着冻融循环次数的增大，拱顶处的体积含水率明显增大，拱侧处的体积含水率略有升高。经过20年的冻融循环，拱顶处与拱侧处的滞水阶段体积含水率分别升高了10.46%与4.21%。（4）随着冻融循环次数的增大，拱脚处与拱底处的体积含水率均略有降低。此为重力场驱动下水分向下渗流，且拱脚处与拱底处无衬砌结构阻隔水分所致。

图8 不同位置处的含水率变化Fig. 8 Water content changes at different positions

图9 为第20 年最大冻深时的温度场分布与衬砌米塞斯应力分布。可以看出：（1）冻结区域在隧道外围呈现差异分布的特征，说明在冻融圈层整体冻胀模型中应考虑冻结厚度分布的不均匀性。（2）顶拱的弧度较大，围岩冻结深度整体较为均匀，约为103 cm。（3）仰拱冻结深度的最大值在其中心位置处，为136 cm。（4）拱脚处的围岩冻结深度最小，仅为88 cm，与最大围岩冻深相差48 cm。此为衬砌几何，以及含水率分布的差异综合导致。一方面拱脚处较大的几何曲率导致相同的冷量要辐射至更大的空间，进而导致该区域冻结深度的减小；另一方面，拱脚较仰拱中部的围岩含水率更高，水-冰相变对冻结锋面的阻滞作用更强。（5）在围岩冻胀力作用下，顶拱的应力分布较为均匀，而仰拱半径较大，在冻胀力作用下发生了一定程度的反向弯曲。（6）结合图9（a）的温度场分布结果，顶拱与仰拱的围岩冻结范围均较拱脚大，使拱脚处发生弯折，最终导致拱脚处衬砌内表面von Mises应力值最大。

图9 最大冻深时的温度场分布与衬砌应力分布Fig. 9 Temperature field and lining stress distributions when the frozen depth reaches largest：temperature distribution （a）， von Mises stress distribution of lining （b）

图10 为第20 年最大冻深时衬砌外围的法向与切向应力分布。其中法向应力与切向应力变量的提取如式（11）所示。可以看出：（1）东天山隧道的衬砌外围法向应力整体上为顶拱＞仰拱＞拱脚。（2）顶拱的弧度整体较大，因此其衬砌外围的法向应力分布较为均匀，同时顶拱的最大法向应力发生在顶拱至拱脚的过渡段。（3）拱脚段衬砌外围的法向应力整体为中间小两端大，局部出现拉应力，最大拉应力为0.031 MPa。（4）仰拱段衬砌外围的法向应力整体为两端大，中间小。（5）除法向应力外，衬砌外围同时存在切向应力。在顶拱、仰拱段，衬砌外围的切向主要承受压应力，其中顶拱两侧、仰拱底部的衬砌外围切向应力分别达到0.228 MPa 与0.462 MPa。（6）拱脚段衬砌外围在切向上承受着一定程度的拉应力。考虑到混凝土材料抗拉强度远小于抗压强度的特性，隧道冻胀力设计中应对该区域予以一定程度的重视，可采用局部增强配筋或布置保温板等方式来保护拱脚段混凝土结构。

图10 最大冻深时的衬砌外围应力分布Fig. 10 Distributions of tunnel frost heave force when the frozen depth reaches largest：normal stress distribution （a）， tangential stress distribution （b）

式中：α为衬砌外法线与x轴的夹角；σN、στ分别为衬砌外围的法向应力与切向应力（MPa）。

对比图10 与图3，可发现青沙山隧道与东天山隧道在拱脚段的法向应力在衬砌外围中的整体分布差异较大。此为水-热-力多场耦合的结果：（1）图9（a）中的温度场可看出，拱脚处的冻结深度最低，若不考虑衬砌几何对围岩冻胀力分布的影响，则拱脚处的法向应力应为最小。（2）实际上，从图9（b）可看出，围岩冻胀力导致衬砌发生了一定程度的变形，顶拱、仰拱的变形导致拱脚处折叠并向外侧围岩挤压，这一效应会增加拱脚段衬砌外围的法向应力。衬砌外围法向应力的分布实际为效应（1）与效应（2）的博弈结果。对于青沙山隧道而言，衬砌厚度较小，效应（2）较效应（1）占优势，因此拱脚段衬砌外围法向压应力较大；而东天山隧道衬砌厚度较大，效应（1）占优势，导致东天山隧道拱脚段的法向压应力较小。

3.3 不同因素对冻胀力分布的影响

3.3.1 最低气温的影响

从隧道洞口至一定纵深处，隧道内部的气温振幅逐渐降低。到达一定隧道纵深时，衬砌背后的围岩基本处于正温范围。基于式（6）中哈密地区的气温变化函数，通过改变气温振幅（其余参数取值与3.2 中相同），计算最低温度分别为-10 ℃、-15 ℃、-20 ℃、-25 ℃下的隧道冻胀力分布。

图11 为第20 年不同最低气温下的隧道冻胀力在极坐标系下的分布，其中顶拱圆心为几何坐标原点。可以看出：（1）衬砌外围的应力分布整体体现为受压。但受温度场差异分布的影响，东天山隧道衬砌外围的法向应力与切向应力均在拱脚处较小甚至局部区域体现为受拉。（2）东天山隧道衬砌外围法向应力整体为顶拱＞仰拱＞拱脚，切向应力整体为仰拱＞顶拱＞拱脚。其中最大法向应力位于顶拱至拱脚过渡段，最大切向应力发生于仰拱中部。（3）随着最低气温的降低，衬砌外围的法向应力与切向应力整体有所增大。其中，顶拱至拱脚过渡段的法向应力、顶拱与仰拱段衬砌外围的切向应力均明显增大。

图11 不同最低气温下的隧道冻胀力分布Fig. 11 Distributions of tunnel frost heave force under various minimum air temperatures：normal stress distribution （a）， tangential stress distribution （b）

3.3.2 初始地层含水率的影响

围岩冻胀为复杂的水-热-力耦合过程，初始地层含水率对于围岩冻胀过程中的热质运移与水-冰相变均存在一定程度的影响。根据东天山隧道实际情况，设置初始地层体积含水率分别为24%、26%、28%、30%（其余参数取值与3.2节中相同），计算相应的冻胀力分布。

图12 为不同初始地层含水率下隧道冻胀力在极坐标系下的分布，其中顶拱圆心为几何坐标原点。可以看出：（1）不同初始地层含水率下，隧道衬砌外围的法向应力与切向应力分布形式相同。其中，顶拱、仰拱段衬砌外围的法向与切向均体现为压应力，拱脚处压应力较小甚至局部轻微受拉。（2）随着初始含水率的增加，顶拱、仰拱段的衬砌外围法向应力略有增大，而顶拱至拱脚过渡段衬砌外围的法向应力明显增大。（3）随着初始地层含水率的增大，衬砌外围切向应力略有增大。（4）在高含水地层，应特别关注顶拱至拱脚过渡段较大的法向压应力，以及拱脚处的局部受拉区域，采取保温或结构增强的措施以减少隧道的冻害。

图12 不同初始地层含水率下的冻胀力分布Fig. 12 Distributions of frost heave force under various initial soil water contents：normal stress distribution （a）， tangential stress distribution （b）

3.3.3 已冻与未冻围岩弹性模量比的影响

图13 为第20 年时不同已冻与未冻围岩模量比下的隧道冻胀力分布。保持未冻结围岩的弹性模量不变，冻结围岩弹性模量分别取未冻结围岩模量的1.2、1.8、2.4、3.0 倍，其余参数与初边值条件与3.2 节中相同。可以看出：（1）不同冻结与未冻围岩模量比下的衬砌外围冻胀力分布规律相同，均在顶拱与仰拱段表现为受压，拱脚段压应力较小且局部受拉。（2）随着冻结与未冻围岩模量比的增大，衬砌外围的法向压应力逐渐减小，这与黄继辉等［20］的研究结果趋势相同。其中顶拱、仰拱处的法向压应力减小幅度最为明显。在模量比由1.2 增大至3.0时，顶拱中部、仰拱中部衬砌外围的法向压应力分别增大了0.03 MPa 与0.01 MPa。（3）冻结与未冻结围岩模量比的增大导致了衬砌外围切向压应力的减小。其中，顶拱与仰拱衬砌外围的切向压应力变化最为明显。在冻结与未冻结围岩模量比由1.2增大至3.0的过程中，顶拱中部、仰拱中部衬砌外围的切向压应力分别降低了0.10 MPa与0.18 MPa。

图13 不同冻土与未冻土模量比下的冻胀力分布Fig. 13 Distributions of frost heave force under various modulus ratios of frozen and unfrozen rock：normal stress distribution （a）， tangential stress distribution （b）

3.3.4 冻胀变形各向异性系数的影响

在常规的室内冻胀试验中，一般将试件置放于完全侧限的容器中进行单向冻结，此时的试件仅在热流方向上发生冻胀变形。而实际工况下，岩土体的冻胀变形主要发生于热流方向，与热流正交方向亦会发生相应的冻胀变形，体冻胀率在空间上分配至热流方向与垂直于热流的方向。因此，围岩冻胀各向异性系数的引入可使得隧道冻胀力的计算与分析结果更加明确。

图14 为不同冻胀变形各向异性系数下的隧道冻胀力分布。可以看出：（1）不同冻胀变形各向异性系数下的衬砌外围冻胀力分布模式相同。其中衬砌外围法向压应力为顶拱＞仰拱＞拱脚。衬砌外围的最大法向压应力发生于顶拱至拱脚过渡段。（2）随着冻胀变形各向异性系数的增大，冻胀应变的方向逐渐集中化，这导致了衬砌外围法向冻胀力的增大。其中顶拱中部衬砌外围的法向压应力增大了0.04 MPa，仰拱中部衬砌外围的法向压应力增大了0.01 MPa。（3）冻胀变形各向异性系数的增大导致仰拱的切向压应力显著提升。在围岩冻胀变形各向异性系数由0.7 增大至1.0 的过程中，仰拱中部衬砌外围的切向压应力由0.46 MPa 增长至0.69 MPa。这是由于冻胀变形方向集中化促进了仰拱的向上弯曲。

图14 不同冻胀变形各向异性系数下的冻胀力分布Fig. 14 Distributions of frost heave force under various orthotropic coefficients of frost heave deformation：normal stress distribution （a）， tangential stress distribution （b）

4 结论

为研究季节冻土区隧道结构的冻胀力分布特征。综合考虑围岩冻胀过程中水分场、温度场与应力场之间的耦合效应，构建并验证了正冻围岩的正交各向异性冻胀模型。并针对东天山隧道，构建了相应的水-热-力耦合数值模型。通过对计算结果分析，得到如下结论：

（1）在大气负温的影响下，衬砌后围岩的温度场、水分场分布并不均匀。由于拱脚处几何曲率最大，相同的冷量需辐射至更大的空间，东天山隧道的围岩冻结深度在拱脚处最小，在仰拱中心处最大，隧道整体的围岩冻结深度差值高达48 cm。单个冻结周期内，围岩的体积含水率变化可分为冻结、融化、滞水、滞水消散4 个阶段。随着冻融循环次数的增加，拱顶处的含水率明显增大，拱侧处的含水率略有升高，拱脚处与仰拱处的含水率则逐年减小。

（2）冻结围岩对衬砌结构同时产生了法向应力与切向应力。东天山隧道法向应力在顶拱两侧与仰拱两端最大，在拱脚处最小且局部体现为拉应力。切向应力在顶拱两侧与仰拱中部最大，在拱脚处最小且局部体现为拉应力。鉴于混凝土的抗拉强度远低于抗压强度，建议在东天山隧道结构设计中，对拱脚处采用保温或结构增强的措施以提升该处的抗冻耐久性。

（3）衬砌外围冻胀力的分布主要受到温度场分布不均匀性与隧道几何的综合影响。拱脚处较小的冻结深度效应导致该处的冻胀力较小，甚至使拱脚处局部产生拉应力。但对于马蹄形隧道而言，衬砌几何使得围岩冻胀力下的顶拱与仰拱发生折叠弯曲，促使拱脚向外侧围岩挤压，这将提高拱脚处的压应力。两种效应的博弈作用决定着隧道冻胀力的分布。

（4）不同最低气温、初始地层含水率、已冻与未

冻围岩模量比、冻胀变形各向异性系数下的隧道冻胀力分布模式整体相同。其中最低气温的降低、冻胀变形各向异性系数的增大分别导致了冻结范围的增大与冻胀应变方向的集中化，明显增大了衬砌外围的冻胀力。已冻与未冻结围岩模量比与隧道冻胀力负相关，而初始地层含水率与隧道冻胀力正相关。