穿黃隧洞縱向變形三維有限元分析

謝小玲，蘇海東，林紹忠（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。