穿黄隧洞纵向变形三维有限元分析

谢小玲，苏海东，林绍忠（1．水利部水工程安全与病害防治工程技术研究中心，武汉 430010；2．长江科学院材料与结构研究所，武汉 430010；3．长江科学院非连续变形分析实验室，武汉 430010）

穿黄隧洞纵向变形三维有限元分析

谢小玲1，2，苏海东1，2，林绍忠1，3
（1．水利部水工程安全与病害防治工程技术研究中心，武汉 430010；2．长江科学院材料与结构研究所，武汉 430010；3．长江科学院非连续变形分析实验室，武汉 430010）

穿黄隧洞深埋于河床摆幅大、冲淤变化剧烈的黄河河道，地处7度地震区，是一个由数百个环拼装联结的非连续变形柔性结构。河床冲淤变化以及地震荷载对隧洞纵向变形以及纵向变形缝的影响问题，是隧洞方案是否可行的关键技术问题之一。为此，采用自编的考虑接触的三维非线性有限元分析程序，计算隧洞结构纵向变形、变形缝张开度、错动位移等。针对隧洞模型计算范围大，接触面多等特点，通过采取多项措施控制计算规模，保证计算精度，结合程序的方程组快速求解方法，解决了穿黄隧洞纵向变形这一技术难题。计算结果表明，在最不利荷载组合作用下，隧洞变形缝的张开度、错动位移和相应跨缝螺栓应力以及混凝土应力均在设计控制的范围内，满足设计要求。

穿黄隧洞；纵向变形；接触非线性；三维有限元

1 概 况

南水北调中线工程穿黄隧洞穿越的河段属典型游荡性河道，河势演变复杂，河床冲淤变化剧烈，主槽摆幅较大。河床主槽摆动产生的冲淤荷载以及地震荷载等因素，可能使隧洞纵向发生不均匀的沉降或回弹变形。而过大的不均匀变形会对衬砌结构、纵向变形缝止水等产生不利影响，因此隧洞设计必须将隧洞纵向的不均匀变形控制在允许的范围内。

对于隧洞这种非连续变形柔性结构的受力特性分析，国内外学者也进行了很多研究，提出了多种计算模型［1，2］，如梁－弹簧模型、梁－接头模型、纵向等效连续化模型以及三维骨架模型等，但这些模型应用于实际工程还存在很多不足，或不能准确模拟管片接头的非线性特性，或模型参数需依据试验选取，或没能反映地基的影响等。文献［2，3］记载，有限元计算中由于隧道纵向长度较大，计算量庞大，造成纵向接头难以模拟，而将隧洞材料均质化，赋予其等效的纵向刚度，使得有限元模型与隧道实际纵向连接情况相差甚远，无法反映隧洞的非连续变形特性。

长江科学院采用自主开发的三维非线性有限元分析程序，计算隧洞结构应力、变形、环缝张开度、缝面接触状态以及螺栓应力等，为隧洞设计提供依据。程序中的接触算法和改进的对称逐步超松弛预处理共轭梯度迭代（SSOR－PCG）的方程组解法，解决了计算模型规模大、接触面多等造成的计算量庞大的关键技术问题。

2 隧洞介绍

穿黄隧洞设计为双线布置，单线长数千米，轴线间距30多米，埋深在30 m以下，穿越的河床地基由亚粘土、粉细砂、中砂、粘土岩组成，其地质系为Q24，Q14，Q2，N，其中亚粘土、粉细砂的力学性质较差。

隧洞衬砌结构设计为双层，分外衬和内衬，外衬由管片拼装成整环，内衬为现浇的整环结构。盾构施工过程中，首先由盾壳撑住已开挖的周围土体，管片在未承载情况下用螺栓将其连接成整环，此后盾构机向前推进，构造下一个环，环与环之间通过螺栓连接形成纵向长数公里的整体隧洞结构。

纵向数个环为一个衬砌段，段与段之间设置纵向变形缝，环向8个管片之间设有环向变形缝，纵、环向变形缝处均设有软垫层，以满足具有一定刚度的柔性结构的需要，衬砌段纵向间距由进口段的1．2 m过渡到主槽范围内标准段的近10 m。

3 计算模拟方法

3．1 计算模型概化

综合河床主槽摆动冲淤荷载的作用范围，计算模型取自进口段竖井边墙，模拟长度约2 500 m。土层模拟范围，垂直向下模拟35 m，向上至地表，自上而下的4个土层，N按实际分布（竖向、纵向）进行模拟。由于垫层（变形缝处）的尺寸单位是mm级，而整个结构模型长2 000多米，两者尺寸悬殊如此大，如何做到保证计算精度，控制计算规模，成为首先要解决的问题。

合适的计算模型和合理的网格布局是问题的关键。分析认为，双洞结构对称，所受的冲淤荷载对称，因此计算模型可只取其中的一洞。试算发现，单个隧洞的变形基本上也呈左右对称分布，因此，仅取单个隧洞的一半进行计算，整体规模减小了3／4。

隧洞纵向的数百个纵向变形缝均需按接触考虑，若再模拟管片间环向变形缝的接触，其计算规模的巨大和多重接触的复杂性将使计算难以进行。考虑到我们研究的重点是纵向变形，因此仅模拟纵向变形缝的接触，并假定混凝土内外衬是一个整体，同时混凝土和土体之间按完好粘结考虑。纵向变形缝处的软垫层（厚5．4 mm）按传压不传拉模拟，垫层和混凝土之间的摩擦系数取0．3，变形缝处沿环向均匀布置的29根32的跨缝螺栓采用杆单元模拟。

隧洞的纵向变形主要反映在变形缝软垫层上，因此，每一段衬砌先沿纵向划分4个单元进行粗算，然后在变形较大的局部区域将网格细化重新计算，由此得到较精确的应力和变形，这样既控制了计算规模，又保证了计算精度。

图1 模型横截面及三维局部网格（三维局部网格沿纵向较粗的线条即缝处）Fig．1 The cross section and 3D local grids in model

模型网格见图1，土体和衬砌均采用8节点的非协调元。模型底部全约束，除顶部表面外，其它边界均施加法向约束，但进口段竖井边墙表面未加约束。

3．2 接触模拟方法

在纵向变形缝处的软垫层采用厚度很小的8节点接触单元［4］对接触面进行模拟，认为缝面能传压，不传拉。

设缝面摩擦系数为f，初始法向间隙为w0。在荷载作用下产生的缝面两侧法向（n）、切向（t，s）的相对位移分别为wr，ur，vr，则缝面接触应力与相对位移之间的关系：

式中：kn，kt，ks分别为缝面单位面积的法向刚度和切向刚度，取值分别为垫层的弹性模量、剪切模量与厚度之比；σn，τt，τs分别为缝面的法向应力和切向应力。

wr＋w0≤0表示法向闭合。如果初始间隙w0＝0，且wr＞0，表示缝面法向张开。当缝面法向张开时，缝面不传递任何应力；当缝面法向闭合时，切向应力可能超过抗剪强度而产生滑移，因此切向应力还要满足条件（2）。

接触计算是非线性问题，一般需要将荷载细分，多步施加，以上一步的缝面接触状态和接触应力作为本计算步的初始值，用变刚度法进行接触问题非线性迭代，直至前后2次的计算结果接近为止，然后转入下一计算步。

3．3 方程组快速求解

模型单元总数为21．2万，节点总数为25万，自由度总数达到75．2万，还要考虑接触非线性迭代，计算量非常庞大。

计算程序中采用对称逐步超松弛预处理共轭梯度迭代法（SSORPCG）的改进迭代格式［5］。与常规的大型有限元方程组的一维变带宽存储的三角形分解直接解法相比，在存储量和计算量方面都降低一个数量级以上。

4 计算参数

4．1 材料参数

材料参数见表1。其中土层的材料参数取自长江科学院的试验结果，按偏于安全考虑，弹模取试验值中的低值。土层弹模为压缩模量，在冲刷荷载段（从荷载图中可见的回弹荷载段）采用卸荷模量，其值为压缩模量的2倍。

表1 材料的力学参数Table 1 Thematerialmechanical parameters

4．2 计算荷载

4．2．1 冲淤荷载

按整体最不利情况而概化的河床主槽摆动冲淤情况见图2。从进口段起始，淤积荷载深度最大约0．9 m，其后冲刷荷载最大冲刷深度约1．3 m。

图2 河床主槽摆动冲淤示意图Fig．2 The riverbed siltingerosion sketch

4．2．2 内水压力荷载

按50 m水头的内水压力计算，计及水重的影响。

4．2．3 地震荷载

隧洞位置相应50年超越概率5%的地震动加速度为a＝0．115 g，采用地震系数法计算，每延米隧洞衬砌结构和内部水体的竖向惯性力幅值为

式中：mq为衬砌结构每延米的质量；mw为洞内水体每延米的质量；g为重力加速度；rH，rB分别为隧洞外半径和内半径。

隧洞所在土层为砂壤土和砂土，地震波在土层中的传播速度为250 m／s，按地震波卓越频率1．33 Hz计算，波长为188 m。假设作用在隧洞上的竖向惯性力沿隧洞纵向按余弦变化，并自竖井侧沿隧洞纵向分布，则2种最不利情况为：

式中z为惯性荷载作用点距竖井侧的距离。

5 成果分析

隧洞纵向变形见图3、图4。

图3 河床冲淤荷载作用下隧洞变形放大示意图Fig．3 The enlarged sketch of tunnel deformation w ith the siltingerosion loading

图4 内水压及地震荷载作用下隧洞变形示意图Fig．4 The tunnel deformation under the action of internalwater pressure and the earthquake loading

冲淤荷载作用下，隧洞在淤积荷载作用段沉降，冲刷荷载段回弹，与该荷载沿纵向的分布形态是相符的。内水压力作用下，隧洞整体下沉，图4中内水压力加地震荷载作用下的隧洞变形，既反映了水重引起的均匀沉降，也反映了地震荷载沿纵向的余弦变化规律。

将相位差为π的2种地震荷载的最不利部位（余弦变化的波峰与波谷）安排在河床摆动冲淤荷载作用下隧洞变形缝最大处，使之同时作用而产生最不利影响，同时考虑内水压力，计算隧洞在冲淤荷载、水压荷载及地震荷载共同作用下的变形，计算结果见表2。

表2 隧洞最大变形Table 2 Themaximum tunnel deformation

由表2可知，冲淤荷载、水压荷载和地震荷载共同作用下的隧洞最大沉降（总变形量）都在180 mm左右，而冲淤荷载引起的隧洞沉降、回弹变形最大值为150 mm和76 mm，沉降变形占总变形量的80%以上，水压荷载引起的沉降量在25 mm左右，占总变形量的15%以下，可见，冲淤荷载是引起隧洞变形的主要荷载。

表3给出冲淤荷载、内水压力荷载及地震荷载作用下的变形缝张开度、错动位移、缝面剪应力以及螺栓应力最大值。图5为第一条变形缝的接触状态和缝面竖向剪应力分布。

表3 变形缝变形、缝面剪应力以及螺栓应力Table 3 The crack deformation，shear stress and bolt stress

图5 冲淤荷载作用第一条变形缝接触状态和剪应力Fig．5 The contact situation and shear stress of the first crack under the siltingerosion loading

计算结果显示，冲淤荷载作用下，淤积和冲刷过渡段，是隧洞由沉降到回弹变形变化最剧烈的区段，也是变形缝张开最大、跨缝螺栓应力最大的部位（见图3中A点），该处变形缝张开约2．7 mm，螺栓拉应力约280 MPa，考虑内水压及地震荷载后，变形缝张开最大值为3．1 mm，小于设计控制值4 mm，相应螺栓拉应力最大值为325 MPa。

变形缝错动位移和剪应力的较大值出现在沉降起始段和由沉降到回弹的平直段，冲淤荷载作用下缝面错动位移约2．2 mm，考虑内水压力和地震荷载后，错动位移最大值增加到约2．8 mm，小于控制值4 mm，满足设计要求。

表4数据显示，冲淤荷载作用下，混凝土衬砌纵向应力最大值为0．55 MPa（拉）、4 MPa（压），即使考虑内水压力和地震荷载的影响，混凝土纵向应力也在设计控制的范围之内。

表4 隧洞混凝土纵向应力最大值Table 4 Themaximum longitudinal stress values MPa

计算结果表明，河床冲淤荷载是影响隧洞纵向变形和应力的主要荷载，在最不利荷载组合（冲淤、内水压力和地震荷载）作用下，隧洞变形缝最大张开度、垂直方向最大错动位移和相应跨缝螺栓最大应力以及混凝土衬砌纵向最大拉、压应力均在设计控制的范围内，满足设计要求。

6 结 语

基于弹性理论的有限元分析已经非常成熟，但对于长达数千米的隧洞纵向变形分析，国内外还没有见到有关研究工作的报道。

本文所做工作，从计算模型的选取到数百个变形缝的接触模拟等等，更贴近实际工程。通过分析变形缝的张开度、错动滑移以及缝面剪应力、混凝土纵向应力，探明了影响隧洞纵向变形的主要影响因素，分析了各因素对隧洞变形的影响范围，给出了结构变形、应力以及变形缝张开度等量值，其分析结果更能为工程所用，解决了穿黄隧洞纵向变形这一技术难题。

本文研究重点是隧洞的纵向变形及纵向应力，若需分析衬砌的环向应力和变形，可选择纵缝张开最不利的隧洞段，采用子模型方法，在模拟纵缝的基础上，详细模拟管片环缝的接触状态和力学行为。

参考文献：

［1］ 张铁柱．盾构隧道纵向计算模型及其问题分析［J］．山西建筑，2008，34（20）：314－315．（ZHANG Tiezhu．Vertical calculation model of shield tunnel ant its problem analysis［J］．Shanxi Architecture，2008，34（20）：314－315．（in Chinese））

［2］ 黄宏伟，臧小龙．盾构隧道纵向变形性态研究分析［J］．地下空间，2002，22（3）：244－251．（HUANG Hongwei，ZANG Xiaolong．Research and analysis on longitudinal deformation characteristics of Shield Tunnel［J］．Underground Space，2002，22（3）：244－251．（in Chinese））

［3］ 徐 凌，黄宏伟，罗富荣．软土地层盾构隧道纵向沉降研究进展［J］城市轨道交通研究，2007，（6）：53－56．（XU Ling，HUANG Hongwei，LUO Furong．Research on longitudinal settlement of shield tunnel in soft ground［J］．Transportation Research on Urban Railway，2007，（6）：53－56．（in Chinese））

［4］ 徐跃之，黄作森，林绍忠．三峡永久船闸第三闸首接触问题仿真计算［J］．长江科学院院报，1998，15（4）：23－27．（XU Yuezhi，HUANG Zhuosen，LIN Shaozhong．Simulation calculation on contact problem fro third head of TGP＇s permanent lock［J］．Journal of Yangtze River Scientific Research Institute，1998，15（4）：23－27．（in Chinese））

［5］ 林绍忠．用预处理共轭梯度法求解有限元方程组及程序设计［J］河海大学学报1998，26（3）：112－115．（LIN Shaozhong．Application of preconditioned conjugated gradientmethod to finite element equations and programme design［J］．Journal of Hohai University，1998，26（3）：112－115．（in Chinese） ）

（编辑：罗玉兰）

Analysis on Longitudinal Deformation of Tunnel Crossed Yellow River by 3D Finite Element

XIE Xiaoling，SU Haidong，LIN Shaozhong
（1．Research Center on Water Engineering Safety and Disaster Prevention of Minist ry ofWater Resources，Wuhan 430010，China；2．Dept．of Material and Structure，Yangtze River Scientific Research Institute，Wuhan 430010，China；3．DDA Center of Yangtze River Scientific Research Institute，Wuhan 430010，China）

The Yellow River tunnel，located at 7 degree seismic region，is deeply buried in the reach where the swing of the riverbed is larger and siltationerosion variation is severe．The tunnel is a discontinuously flexible structure composed of several hundreds rings．The key technique problem in the feasible scheme is the longitudinal deformation caused by erosionsiltation change and earthquake load．So the 3D nonlinear finite element program is developed to analyze the longitudinal deformation，crack opening，dislocation etc．．According to the characteristics of the large scale model and many contact interfaces，several approaches are adopted to control calculation scale and to ensure calculation precision，including rapid solution algorithm for the equations．The results show that the crack opening，the dislocation displacement，the bolt stress of span joints，and the concrete stress are all in the design permission range，and can satisfy the design requirement．

Yellow River Tunnel；longitudinal deformation；contact nonlinearity；3D finite element

TU311．4；TV314

A

1001－5485（2010）07－0060－05

20091209；

20100407

“十一五”国家科技支撑计划项目“复杂地质条件下穿黄隧洞工程关键技术研究”（2006BAB04A11）

谢小玲（1959），女，广东潮汕人，高级工程师，主要从事水工结构研究，（电话）02782829754（电子信箱）xiexiaoling01＠163．com。