復雜隧洞襯砌結構的三維有限元分析及結構配筋

董正中，肖 明，楊 瑩(武漢大學水資源與水電工程科學國家重點實驗室，武漢 430072)

0 引 言

隨著我國大型水電站建設的展開，水電站隧洞襯砌結構向著體形大型化，地質條件多樣化，受力條件復雜化的方向發展[1]。襯砌自身的結構形式也日益復雜。如何分析襯砌結構的力學特性，以及如何計算襯砌結構的配筋，使之同時滿足安全性和經濟性的要求，是一個重要的課題。壓力隧洞的傳統計算方法[2,3]假定圍巖及襯砌材料為連續各向同性的彈性介質，忽視了圍巖及襯砌材料的非線性性質的不利影響，并且是將圍巖與襯砌結構分開進行計算，從而導致隧洞混凝土襯砌厚度及配筋量偏大。因此，有必要對考慮隧洞圍巖、襯砌聯合承載及材料非線性情形下的襯砌結構進行分析計算。基于此，本文首先采用三維彈塑性損傷有限元分析方法，對復雜受力條件下襯砌結構進行計算，研究探討了各工況下隧洞襯砌結構的受力特點；其次利用有限元插值應力場在襯砌截面上積分計算內力，基于內力結果在配筋截面上進行配筋計算，以保證襯砌配筋滿足要求，并提出優化建議。

該方法基于三維有限元計算，綜合考慮了初始地應力場，開挖卸荷，復雜地質條件及地下水位等多種因素，研究了襯砌結構的穩定特性，并直接由有限元結果插值進行配筋計算。為在工程實際中選擇經濟合理的襯砌形式以及提高隧洞運行的安全度提供依據。相較于傳統襯砌結構設計評價方法更高效，精確，可為類似工程提供借鑒與參考。

1 襯砌結構三維有限元分析的基本方法

本文對襯砌結構采用彈性模型進行計算；對圍巖單元則采用彈塑性模型[4]，屈服準則采用Ziekiewicz-Pande準則，以增量變塑性剛度法進行迭代計算[5]。

襯砌結構是在隧洞開挖完畢后施加上去的 。為了反映隧洞開挖變形對襯砌作用的影響，可以洞室開挖釋放的荷載:

(1)

改寫成：

(2)

式中：σ0是開挖單元的初始地應力場，包括自重應力場和構造應力場；γ為巖體容重；α為荷載分配系數。

襯砌結構施加前，作用于結構的荷載為

(3)

襯砌結構施加后，作用于結構的荷載為：

(4)

式中荷載分配系數α的取值與巖性、地應力和支護時機等因素相關，工程中一般靠經驗確定[6]。因此，本文嘗試給出一種確定取值的數值計算方法。

彈塑性有限元計算中，開挖后圍巖單元的應力狀態按下式可分為彈性和塑性，即：

(5)

式中：F為屈服函數；σ為開挖后的圍巖應力；σ0為巖體的初始應力；Δσ為開挖引起的應力增量。當巖體開挖后，若圍巖單元進入塑性狀態，則一定存在一個臨界應力狀態，滿足：

F({σ})=F({σ0}+β{Δσ})=0

(6)

式中：β為單元的彈性系數。令p=1-β，p稱為塑性系數，其大小反映了總應力增量Δσ中塑性荷載的比例。由式(6)可以看出，初始地應力場下已屈服的巖體單元p=1；開挖后仍處于彈性的單元p=0。

根據式(6)，采用Zienkiewicz-Pande屈服準則[7]，通過一次開挖計算，令全部開挖釋放荷載作用于圍巖，求出所有圍巖單元的塑性系數p；塑性系數p的分布在一定程度上反映了圍巖的承載條件，因此，式(2)中的荷載分配系數 可根據洞室頂拱和邊墻單元的塑性系數p分布按式(7)綜合確定。

α=(1-η)p

(7)

其中，η為支護時機滯后系數。視工程實際情況，對于自承能力較差的巖體，數值分析時可取用較小η值，使圍巖單獨承擔的荷載減少；對于自承能力較強的圍巖，數值分析時可適當放大η取值，使圍巖承擔部分塑性荷載。一般地，巖性越差，地應力越大，則η越小；反之，η越大。

洞室開挖荷載釋放完畢后，再根據運行期引水隧洞的水位分布情況，計算隧洞襯砌受內、外水壓力作用的受力特征 。

2 襯砌結構的配筋計算

本文基于應力圖形法，編制了地下洞室襯砌結構的配筋計算程序。并可以自動生成配筋截面：對于襯砌結構有限元模型，先生成模型的拓撲關系，即結點、單元線、單元面、單元間的相互包含關系；再對模型進行消隱，建立表面單元線表，這里定義表面單元面僅有一相關單元或相關單元僅有一襯砌單元，表面單元線包含在表面單元面中，以一平面截取表面單元線，生成內外輪廓交點，依序連接交點，即形成內外輪廓線，即可獲得配筋截面。

2.1 應力修正與內力修正

由于有限元計算的特點，單元內部有應力均化現象，從而導致配筋截面上計算彎矩偏小。引入應力修正系數r。首先對結點應力作估值，可由與該結點相關的單元應力取平均值求得。分別計算單元應力及結點應力在問題域局部的極值，令其為σeij、σnij，取

r=‖σeij‖/‖σnij‖

(8)

這里‖ ‖表示應力張量的度量，取為應力不變量的函數。一般有r>1。將r與結點應力相乘即得修正結點應力。

有時沿配筋截面僅有少數幾層單元，這可能導致計算彎矩偏小甚至方向相反。引入應變修正方法，假設在配筋截面方向上，應力依線性分布，設應力梯度為沿配筋截面應力的斜率，則應力梯度可由配筋截面兩端點應變計算獲得。假設配筋截面兩端點沿所取投影方向正應變分別為ε1和ε2，則可得配筋截面上應力梯度為：

k=γ(ε1-ε2)/L

(9)

式中：L為配筋截面長度；L依結構形狀和受力情況取值。

2.2 配筋截面的內力及配筋計算

單個截面的內力計算過程如下：

(1)設置應力插值點，取配筋截面分點作為插值點，份數與單元層數成比例，一般為其2倍或以上[8]。

(2)由上述插值方法計算各插值點應力，并根據上述應力修正方法對應力進行修正。

(3) 由柯西公式σn=ninjσij計算各插值點沿某一方向正應力。通過數值積分計算配筋截面上內力，這里采用復合積分公式。依據沿配筋截面方向單元層數對截面上內力作修正[9,10]。

以一長圓筒為例。其內徑r=5 m，外徑R=6 m，受外壓q1=0.5 MPa，內壓q2=0.7 MPa，按彈性本構計算，則由解析解公式可得配筋截面上軸力N=500 kN，彎矩N=16.56 kN·m。建立該圓筒有限元模型，沿徑向取5層單元，對其進行有限元分析。圖1左右分別為經由插值及數值積分所得的軸力及彎矩結果。可以看出，與解析解相比，軸力誤差在4%以內，彎矩誤差在6%以內。若考慮應力修正取應力張量度量為第一應力不變量可得應力修正系數為1.01，從而使結果更趨近于解析解。

圖1 圓筒模型截面軸力彎矩圖Fig.1 Axial force and bending moment graph of cylinder model

截面上配筋參照《水工混凝土結構設計規范》(DL/T 5057-2009)。如果截面上應力分布接近線性，按正截面承載力方法計算配筋；若偏離線性較大，按拉應力圖形面積計算，即As=KT/fy，這里T取為拉應力面積， 為鋼筋強度設計值，K為承載力安全系數。裂縫寬度驗算依據具體工程情況參照相應規范計算。

3 工程實例分析

3.1 工程概況

本文結合某水電站引水隧洞襯砌結構進行分析計算。電站裝機容量2 400 MW。電站單機引用流量較大(Q=621.4 m3/s)，機組采用單管單機供水形式，對應4臺機組。每條引水道由進口漸變段(矩形→圓形)、上平段(有壓引水隧洞)、漸縮段(有壓隧洞段→地下壓力鋼管段)、上彎段、斜井段、下彎段、下平段等組成。引水道進口底板高程1 575.00 m，出口中心高程為1 494.80 m。由于引水隧洞條數較多、洞徑較大(最大開挖斷面16.8 m，居已建發電引水隧洞最大開挖直徑的前茅)、相鄰隧洞間的巖柱厚度較小(略大于一倍開挖洞徑)，且圍巖地質條件較復雜，因而有必要對隧洞襯砌結構受力特性進行分析，并進行結構配筋計算以滿足限裂要求。

3.2 計算參數及有限元模型

計算網格一共剖分了362 544個8節點空間等參單元，節點總數376 960個。由四條引水隧洞進口段始端建至有壓段末端。有限元模型參見圖2，根據引水隧洞地址剖面圖，共劃分了三類巖體單元。襯砌采用C25混凝土，襯砌結構參見圖3。計算所取的巖體和襯砌參數見表1及表2。

圖3 襯砌單元模型Fig.3 Lining element model

表1 地下洞室圍巖參數Tab.1 Surrounding rock parameters of cavern

表2 襯砌材料物理力學參數設計值Tab.2 Design values of mechanical parameters of lining

模型X軸垂直于進水口段軸線方向，Y軸沿進水口段水的流向，Z軸與大地坐標重合，指向上為正，范圍由1 433.6 m高程至地表。邊界條件為：底部全約束，四周法向約束，頂部不約束。

3.3 計算工況

本文以四條隧洞均過流為最不利工況展示襯砌受力及配筋成果。

(1)內水壓力由正常蓄水位1 618 m確定。水擊壓力，4號壓力管道末端的最大水擊壓力值為36.8 m，其余3個引水道的壓力管道末端的水擊壓力值也按4號引水道取，為36.8 m，水擊壓力分布按線性規律計算。

(2)外水壓力由水庫蓄水后地下水位線確定，水庫蓄水后的地下水位線由滲流場分析求得[11,12]。根據計算所得的巖體滲流場節點水頭，在襯砌結構上施加滲透體積力。

(3)計算使用自行編制的地下洞室三維有限元分析以及滲流分析程序。

3.4 計算結果分析

(1)襯砌應力及位移分布規律。進水口漸變段1-1以及上彎段2-2兩段截面襯砌第三主應力分布規律參見圖4及圖5。兩段典型截面變形示意圖參見圖6及圖7。

受內水壓力作用，隧洞典型斷面第一主應力均為壓應力，分布在-1.14～-0.19 MPa，應力矢量沿洞周徑向。

圖4 襯砌1-1斷面第三主應力分布圖Fig.4 3rd principal stress of 1-1 section of lining

圖5 襯砌2-2斷面第三主應力分布圖Fig.5 3rd principal stress of 2-2 section of lining

圖6 襯砌1-1斷面變形示意圖(單位：mm)Fig.6 Displacement of 1-1 section of lining

圖7 襯砌2-2斷面變形示意圖(單位：mm)Fig.7 Displacement of 2-2 section of lining

斷面第三主應力為拉應力，應力矢量沿洞周切向。1-1斷面由于斷面形式，四周拐角處出現應力集中，拉應力值達1.5 MPa；2-2斷面由于高程比1-1斷面高程低，內水壓力較大，斷面各部位第三主應力量值均大于1-1斷面第三主應力。4號洞由于洞徑最大，洞周圍巖參數較其他三條隧洞洞周圍巖參數低，故而襯砌承擔較多圍巖荷載，斷面各部位拉應力量值大于其他隧洞拉應力值，兩側腰部拉應力值達2.5 MPa。計算結果表明，漸縮段以及下游襯砌拉應力超出混凝土抗拉強度設計值，應重點加強配筋。

1-1斷面位移分布在0.2～0.7 mm。2-2斷面各部位位移值相對于1-1斷面位移略大，分布在0.3～1.2 mm。變形指向洞外，各斷面底部位移大于隧洞頂部及兩側位移。其中圓形斷面洞周變形較漸變段洞周變形均勻。

(2)隧洞襯砌配筋計算限于篇幅，本文展示4號隧洞1-1斷面和2-2斷面的截面內力以及配筋分布圖，分布參見圖8和圖9。

圖8 4號隧洞1-1斷面配筋示意圖Fig.8 Reinforcement drawing of 1-1 section of lining

圖9 4號隧洞2-2斷面配筋示意圖Fig.9 Reinforcement drawing of 2-2 section of lining

1-1斷面受截面形狀影響，角點處彎矩與邊上彎矩方向相反，由于內水壓力，邊上彎矩向內彎，近角點處彎矩向外彎。其余各斷面彎矩基本向外。

各洞基本按小偏心受拉配筋，并利用裂縫寬度限值驗算配筋率。按所處環境類別，取裂縫寬度限值為0.3 mm，先由單寬配筋截面上內力進行配筋計算，結果表明：1-1斷面按構造配筋，2-2斷面裂縫寬度均未超出限值。若對于其他斷面出現裂縫寬度超出限值，需對這些截面通過裂縫寬度驗算求得配筋率。

4 結 論

(1)考慮隧洞圍巖、襯砌聯合承載及材料非線性情形下對襯砌結構進行分析計算。為了反映隧洞開挖荷載對襯砌作用的影響，提出了“荷載分配系數”取值的一種數值計算方法，為類似工程計算提供參考。

(2)基于有限元分析結果，實現了地下洞室隧洞襯砌結構的配筋。采用應力圖形法配筋，配筋截面自動生成，能適應各種截面形式，同時應力分析與配筋過程相對獨立，具有較高的效率，對圓筒結構的驗算表面，該方法具有較高的精度。

(3)采用應力修勻計算結點應力，再計算插值點應力。特別地，該方法對于網絡稀疏的情況具有很好的精度。對某水電工程引水隧洞的計算實例表明，襯砌結構有限元分析以及配筋計算方法成果合理，可以用于實際工程的配筋設計中。

□

[1] 卞 康, 肖 明, 任 祎. 復雜調壓井結構的三維有限元分析計算[J]. 中國農村水利水電, 2008,(3):115-117.

[2] 河海大學, 武漢大學, 大連理工大學,等. 水工鋼筋混凝土結構學[M]. 北京：中國水利水電出版社, 2009.

[3] 呂西林, 金國芳, 吳曉涵. 鋼筋混凝土結構非線性有限元理論與應用[M]. 上海：同濟大學出版社, 1997.

[4] 肖 明. 地下高壓岔管洞室穩定及襯砌損傷開裂分析[J]. 武漢水利電力大學學報, 1995,28(6):594-599.

[5] 俞裕泰, 肖 明. 大型地下洞室圍巖穩定三維彈塑性有限元分析[J]. 巖石力學與工程學報, 1987,6(1):47-56.

[6] 肖 明. 地下洞室施工開挖三維動態過程數值模擬分析[J]. 巖土工程學報, 2000,22(4):421-425.

[7] 王勖成, 邵 敏. 有限單元法基本原理和數值方法[M]. 北京：清華大學出版社有限公司, 1997.

[8] 張雨霆, 肖 明, 熊兆平. 三維空間離散點數據場的插值方法[J]. 武漢大學學報: 工學版, 2008,41(4):34-37.

[9] 錢向東, 任青文, 趙 引. 一種高效的等參有限元逆變換算法 Ξ[J]. 計算力學學報, 1998,15:4.

[10] 李春光, 鄭 宏, 葛修潤, 等. 六面體單元等參逆變換的一種迭代解法[J]. 巖土力學, 2004,25(7):1-050.

[11] 張 巍, 肖 明. 地下工程復雜滲流場數值模擬研究與應用[J]. 巖土工程技術, 2005,19(1):13-16.

[12] 張 巍, 肖 明. 地下工程滲流斷層數值模擬的隱式復合材料單元法研究[J]. 巖土工程學報, 2005,27(10):1 203-1 206.