基于有限元模擬的跨斷層輸水管道失效風(fēng)險評價
——以阿紅干渠管道為例

顧世祥，梅 偉，唐 暢，楊 帆，霍玉國，曹子君

（1.云南省水利水電勘測設(shè)計研究院，昆明650021；2.武漢大學(xué)水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室工程風(fēng)險與防災(zāi)研究所，武漢430072）

0 引 言

滇中引水工程跨越金沙江、瀾滄江、紅河、南盤江四大水系，穿越橫斷山系高中山地貌區(qū)及滇中、滇東盆地山原區(qū)，沿線涉及多條地震帶?；顒訑嗔褞话l(fā)地震帶來的斷層錯動會導(dǎo)致埋地管道受力發(fā)生“S”型變形，管道可能因發(fā)生局部拉伸破壞或壓縮破壞而失效，不利于整個輸水管道工程的安全運(yùn)行，甚至?xí)韲?yán)重的次生災(zāi)害［1］。研究典型跨斷層埋地管道響應(yīng)能夠?yàn)楸Ｕ匣顒訑嗔褞л斔艿腊踩峁﹨⒖家罁?jù)，對輸水管道工程防震減災(zāi)具有十分重要的意義。

跨斷層埋地管道分析方法包括理論解析法和數(shù)值模擬法。解析方法通常以索模型或梁模型為基礎(chǔ)進(jìn)行研究，該方法難以分析管道受壓和管截面大變形的情況［2］。數(shù)值模擬方法是研究埋地管道地震響應(yīng)的重要手段，管-土相互作用建模方式是數(shù)值模擬方法中的重要環(huán)節(jié)，主要包括兩種方法：土彈簧法和接觸面法。Vazouras 等［3，4］采用非線性接觸面，用殼單元和實(shí)體單元分別模擬管道和土體，對管道在走滑斷層作用下進(jìn)行了精細(xì)化數(shù)值模擬。目前，國內(nèi)外埋地管道抗震設(shè)計指南（例如GB/T50470-2017［5］；ALA，2005［6］；PRCI，2004［7］）推薦采用土彈簧法進(jìn)行跨斷層埋地管道響應(yīng)分析。土彈簧法相比于接觸面法具有計算時間短、收斂性好及應(yīng)用廣泛的特點(diǎn)［8］。ABAQUS 軟件中基于土彈簧法開發(fā)的PSI 接觸單元能夠合理地表征管-土相互作用的非線性行為，在埋地管道和周圍土體的相互作用分析中廣泛應(yīng)用［9-11］。

實(shí)際工程中，斷層錯動引起的管道響應(yīng)受到多種不確定性因素的影響［12］，比如管材力學(xué)參數(shù)的不確定性、土體參數(shù)的不確定性以及荷載不確定性等?；诿商乜宓碾S機(jī)模擬方法能夠合理地考慮上述不確定性對管道結(jié)構(gòu)響應(yīng)的影響。例如，樊恒等［13］研究了考慮地震水平加速度、管壁厚度、屈服強(qiáng)度、彈性模量和泊松比的不確定性對管道動力響應(yīng)的影響。陶凱爾等［14］采用解析法分析跨斷層埋地管道響應(yīng)時考慮了地震烈度和管壁厚度的不確定性，采用數(shù)值模擬分析跨斷層埋地管道響應(yīng)時僅考慮了斷層錯動位移的不確定性。目前，國內(nèi)外學(xué)者對跨斷層埋地管道的不確定性分析研究相對較少，且一般根據(jù)跨斷層埋地管道解析法進(jìn)行分析，主要考慮管材參數(shù)和荷載的不確定性，沒有考慮場地土體參數(shù)的不確定性，研究算例常為走滑斷層。

本文以滇中引水二期工程玉溪段阿紅干渠管道五段為例，研究了活動斷層作用下管道的失效風(fēng)險。利用ABAQUS 軟件中的PSI單元模擬管-土相互作用的非線性行為，通過等效樣本法［15］獲得場地的土體力學(xué)參數(shù)樣本值，考慮了地震位移荷載、管道材料強(qiáng)度以及土體不排水剪切強(qiáng)度的不確定性對管道動力響應(yīng)的影響，運(yùn)用蒙特卡洛模擬方法進(jìn)行隨機(jī)模擬分析，得到所研究管道的失效概率用以評估管道的安全性和可靠性。在此基礎(chǔ)上，研究了地震震級、管壁厚度和填土深度對于管道失效概率的影響規(guī)律，為該工程及類似工程提供參考依據(jù)。

1 場地概況

滇中引水工程二期工程玉溪段輸水管道布設(shè)不可避免地穿越地震活動帶，管道抗震問題突出。斷層活動引起的永久地面變形的發(fā)生場地相對固定，管道損壞率高，管道抗剪斷問題成為管道工程中的重要問題之一。如圖1所示，玉溪段阿紅干渠管道同時跨越普渡河斷裂帶和曲江斷裂帶，其中阿紅干渠管道五段跨越普渡河斷裂帶。

普渡河斷裂帶是滇中“南北向構(gòu)造帶”的主要成分之一?？傮w走向近南北向，傾向向東，傾角約75°～80°［16］，區(qū)域長度大于100 km，斷裂性質(zhì)為逆斷層，左傾滑。其主干斷裂大致沿普渡河延伸，總體走向近南北，全長約300 km。研究表明，普渡河斷裂第四紀(jì)以來的平均水平運(yùn)動速率為0.9～2.0 mm/a，平均垂直運(yùn)動速率為0.35～0.47 mm/a［17］。該斷裂強(qiáng)烈活動的最新時代為早中更新世時期，晚更新世以來，該斷裂帶仍有活動，但活動性明顯減弱。普渡河斷裂帶自公元1600年記錄以來，5 級以上地震17 次，包括1985年4月18日在云南祿勸-尋甸交界地區(qū)發(fā)生的Ms6.3級地震。

根據(jù)玉溪盆地鉆孔探測成果圖［18］，土層深度5.5 m 以上為黏土層，土層深度5.5～13.3 m 為粉質(zhì)黏土。根據(jù)阿紅干渠地質(zhì)資料，黏土層的內(nèi)摩擦角為12°～15°。阿紅干渠管道五段與普渡河斷裂帶相交，夾角約為45°，該管段分配水頭15.298 m，設(shè)計流量為5.9 m3/s，設(shè)計管徑1.8 m，距離管軸的填土埋深為3.2 m，輸水管道擬采用Q235C鋼材。

2 確定性分析模型

2.1 位移分量

在研究跨斷層埋地管道響應(yīng)時，通常將其近似簡化為以位移荷載控制的靜力問題來分析，因此合理地估計斷層位移，對研究特定場地處管道響應(yīng)具有重要作用。本文僅考慮普渡河斷裂帶為逆斷層的情況，根據(jù)地震震級，采用逆斷層地震位移估算經(jīng)驗(yàn)公式［19］，可以得出斷層位移估計值。

式中：M為地震震級；MD為最大地表斷層位移，m；AD為平均地表斷層位移，m；σlog（MD）和σlog（AD）分別為兩個回歸方程的標(biāo)準(zhǔn)差。

需要注意的是，公式（1）和（2）僅適用于平行于斷層走向的豎直平面上的位移，此外其只適用于M為5.4～7.4的情況。在本文有限元模擬中，僅考慮斷層最大位移MD的工況，即：

式中：DFS為平行于斷層走向的豎直平面上的斷層位移，m。

根據(jù)位移與夾角的幾何關(guān)系，可由平行于斷層走向豎直平面上的斷層位移推得斷層的總位移［7］：

式中：Δ為斷層總位移；α為斷層傾角；β為斷層走向與區(qū)域應(yīng)力方位角之間的水平夾角（對于以逆斷層為主的斷層β=90°）。

圖2為典型逆斷層位移分量示意圖。如圖2（a）表示逆斷層總位移與豎直方向位移的關(guān)系，圖2（b）為水平面上平行于管軸方向位移與垂直于管軸方向位移的關(guān)系。已知逆斷層傾角、斷層與管道的夾角，可將斷層總位移分解為三個相互垂直方向的位移［7，20］：

式中：κ為管道與斷層之間的夾角；ΔX為平行于管道軸線方向的位移分量；ΔY為水平面上垂直于管道軸向方向的位移分量；ΔZ為豎直方向的位移分量。

2.2 管-土相互作用

管-土相互作用是影響埋地管道響應(yīng)的關(guān)鍵因素之一。土彈簧法以彈性地基梁假設(shè)為基礎(chǔ)，將土體對管道的作用簡化為3 個方向上的土彈簧，通過土彈簧的特性曲線（圖3）來表征管-土相互作用的本構(gòu)關(guān)系［21，22］。

管-土相互作用的本構(gòu)關(guān)系可以通過室內(nèi)試驗(yàn)獲得，也可以采用文獻(xiàn)中的簡化模型，其中美國ASCE 生命線工程技術(shù)規(guī)程中提出的管－土相互作用模型被廣泛應(yīng)用。本文中管道所埋地層為黏土層，表征土彈簧特性的特征值可計算如下［23］。

式中：D為管道外徑；Su為土體不排水剪切強(qiáng)度；H為距離管道軸線的土體深度；α為與Su有關(guān)的系數(shù)；Nch和Ncv為與H/D有關(guān)的系數(shù)；Nc為與土體內(nèi)摩擦角φ有關(guān)的系數(shù)；tu為管軸方向單位長度土體對管道的最大抗力；xu為yu對應(yīng)的最小軸向位移；pu為水平方向單位長度土體對管道的最大抗力；yu為pu對應(yīng)的最小水平方向位移；qu1和qu2分別為豎直向下和豎直向上方向單位長度土體對管道的最大抗力，zu1為qu1對應(yīng)的最小豎向位移，zu2為qu2對應(yīng)的最小豎向位移。

因此，當(dāng)已知Su、φ、H和D時，即可獲得土彈簧特性曲線。ABAQUS 中內(nèi)置的PSI 單元基于土彈簧模型表征管-土相互作用。當(dāng)采用PSI 單元模擬埋地管道和土體的相互作用時，可采用梁單元、管單元或彎接頭單元來模擬管道，PSI 的一側(cè)與管道共節(jié)點(diǎn)，另一側(cè)代表土體表面，并用于施加邊界條件，以描述地面的運(yùn)動［9，10］。

2.3 建模過程

玉溪段阿紅干渠管道五段管道壁厚為12 mm，外徑為1.824 m，埋深H為3.2 m，φ取14°，初設(shè)土體不排水剪切強(qiáng)度Su為66 kPa。可估算由水壓引起的管道軸向應(yīng)力［24］：

式中：p為管道內(nèi)部水壓；δ為管壁厚度；σn為由內(nèi)部水壓引起的管道軸向拉應(yīng)力。

初步計算可得σn=5.7 MPa，該值遠(yuǎn)遠(yuǎn)小于管材的屈服強(qiáng)度（235 MPa）。因此，建模時可忽略內(nèi)部水壓的影響。采用簡化的折線彈塑性模型表征管材特性，圖4所示為Q235C 鋼材的理想應(yīng)力-應(yīng)變曲線，鋼材的屈服強(qiáng)度為235 MPa，彈性模量為206 GPa，泊松比為0.3，拉伸強(qiáng)度為395 MPa。本文所建立的三維模型如圖5所示，選取管道長度600 m，右側(cè)為斷層上盤，斷層錯動時，右盤上移。建立模型時，管道、地面和管-土相互作用分別采用PIPE31 單元、RB3D2 單元和PSI34 單元進(jìn)行模擬，靠近斷層的位置管單元網(wǎng)格劃分較密，遠(yuǎn)離斷層的位置管單元網(wǎng)格劃分的較稀疏，以便獲得良好的精度和計算效率。根據(jù)2.1節(jié)計算得到地震震級為6.3時的位移荷載分量，如表1所示。

表1 地震工況下的位移分量（M=6.3）Tab.1 Displacement components under earthquake of M=6.3

2.4 失效準(zhǔn)則

基于應(yīng)變的失效準(zhǔn)則相比應(yīng)力失效準(zhǔn)則更能夠發(fā)揮材料的塑性和延展性能，本文基于應(yīng)變設(shè)計方法對通過活動斷層的埋地管道進(jìn)行抗拉伸和抗壓縮驗(yàn)算。管道軸向容許拉伸應(yīng)變和管道軸向容許壓縮應(yīng)變應(yīng)分別按式（13）和式（14）計算［5］：

式中：［εt］F和［εc］F分別為埋地管道抗斷的軸向容許拉伸應(yīng)變和容許壓縮應(yīng)變；εtcrit和εccrit分別為管段的極限拉伸應(yīng)變和壓縮應(yīng)變；φεt為拉伸應(yīng)變承載系數(shù)，當(dāng)環(huán)向應(yīng)力小于屈服強(qiáng)度的40%時，取0.9；φεc為壓縮應(yīng)變承載系數(shù)，取0.6。本文中由GB/T50470-2017［5］附錄D 估算Q235C 設(shè)計容許拉伸應(yīng)變?nèi)≈禐?.91%，設(shè)計容許壓縮應(yīng)變?nèi)≈禐?.55%。

2.5 確定性結(jié)果分析

管道跟隨土體發(fā)生3 個方向的位移錯動，但以豎直方向的位移為主，右側(cè)管道向上運(yùn)動時，管道上部受拉，下側(cè)受壓，由于逆斷層在水平方向的位移分量使得管道軸向受壓，因此該管道可能發(fā)生局部嚴(yán)重受壓。管道軸向應(yīng)力云圖如圖6所示，斷層錯動發(fā)生后，管道發(fā)生“S型”形變，且大部分區(qū)域發(fā)生軸向壓縮應(yīng)變。圖7中紅色區(qū)域表示管道發(fā)生了塑性變形，可見位于斷層附近的管道大部分區(qū)域已進(jìn)入塑性階段。圖8為管道軸向應(yīng)變包絡(luò)圖。如圖所示管道最大拉應(yīng)變εtmax為0.3%，小于［εt］F，管道最大壓應(yīng)變εcmax為2.2%，大于［εc］F，管道因而壓縮失效。

3 不確定性分析

3.1 隨機(jī)變量

跨斷層埋地管道響應(yīng)受到多種不確定性因素的影響。本文中考慮了管材力學(xué)參數(shù)、位移荷載以及土體參數(shù)的不確定性對管道響應(yīng)的影響。鋼材由于其在鍛造過程中存在缺陷或者不均勻性等因素使得材料強(qiáng)度具有不確定性，已知Q235C 鋼材屈服強(qiáng)度235 MPa。本文中假設(shè)鋼材屈服強(qiáng)度服從對數(shù)正態(tài)分布，其均值為235 MPa，變異系數(shù)為10%，即標(biāo)準(zhǔn)差為23.5 MPa，其他強(qiáng)度參數(shù)（如極限抗拉強(qiáng)度與屈服強(qiáng)度的比值）保持不變。地震震級與最大地表斷層位移之間并非滿足映射關(guān)系，其關(guān)系往往由統(tǒng)計分析得到。因此，利用經(jīng)驗(yàn)公式（1）計算最大地表斷層位移時具有不確定性，其中l(wèi)og（MD）服從均值為（-1.84+0.29M），標(biāo)準(zhǔn)差為0.42的正態(tài)分布。

巖土工程勘探數(shù)據(jù)通常十分有限，以至于無法直接根據(jù)這些勘探數(shù)據(jù)獲得有效的土體參數(shù)統(tǒng)計特征和概率分布。基于貝葉斯理論的等效樣本法［15］能夠?qū)龅赜邢薜目碧綌?shù)據(jù)和先驗(yàn)信息結(jié)合，得到所需巖土體參數(shù)的后驗(yàn)分布，并采用馬爾科夫蒙特卡洛模擬方法獲得后驗(yàn)分布的大量等效樣本。

本文采用玉溪盆地地區(qū)黏土層的剪切波速數(shù)據(jù)［18］計算管道所埋地層的不排水剪切強(qiáng)度。不排水剪切強(qiáng)度Su與剪切波速Vs具有以下經(jīng)驗(yàn)關(guān)系［25］：

式中：σs為上述回歸方程中不排水剪切強(qiáng)度Su的標(biāo)準(zhǔn)差。等效樣本法已經(jīng)內(nèi)置于BEST 軟件［15］，本文采用BEST 軟件基于5 個不同深度處剪切波速Vs數(shù)據(jù)獲得該土層30 000 個不排水剪切強(qiáng)度Su的等效樣本，樣本均值為66 kPa，標(biāo)準(zhǔn)差為14.5 kPa。由于BEST 采用馬爾科夫鏈蒙特卡洛模擬產(chǎn)生等效樣本，不同樣本間存在一定相關(guān)性。為了降低樣本相關(guān)性的影響，本文從30 000個樣本中等間距選取1 000個樣本進(jìn)行管道響應(yīng)分析。

在對跨斷層埋地管道響應(yīng)分析進(jìn)行隨機(jī)模擬時，根據(jù)屈服強(qiáng)度和最大地表斷層位移的概率分布生成1 000 組隨機(jī)樣本，這些樣本與Su的1 000 個等效樣本共同作為管道分析的輸入?yún)?shù)，進(jìn)行1 000 次管道響應(yīng)確定性分析，統(tǒng)計管道響應(yīng)分析結(jié)果。

3.2 模型對比

由于材料非線性、大變形和管-土非線性接觸等問題，基于有限元的隨機(jī)模擬中不可避免地存在不收斂的情況。為避免或減少不收斂的情況出現(xiàn)，本文在隨機(jī)模擬中進(jìn)行了模型簡化，由于本算例中基于軸向應(yīng)變失效準(zhǔn)則判斷管道是否失效，軸向應(yīng)變主要與軸向拉壓和管道彎曲有關(guān)。如圖2所示，軸向位移ΔX會直接引起管道軸向拉壓變形，剪切位移ΔY和ΔZ會引起管道彎曲和軸向拉伸。管道在剪切作用下，彎曲變形使得管道兩側(cè)分別承受大小相等的拉壓應(yīng)變，軸向拉伸增大了受拉一側(cè)的拉應(yīng)變，而削弱了受壓一側(cè)的壓應(yīng)變。因此，代表剪切作用的ΔY和ΔZ對拉應(yīng)變的影響較大，而對壓應(yīng)變的影響較小。

由公式（5）～（7）計算可得，軸向位移ΔX和水平向位移ΔY相等，且斷層錯動引起的豎直方向位移ΔZ為軸向位移ΔX的5.3 倍。因此，隨機(jī)模擬中僅考慮施加軸向和豎直方向的位移荷載以簡化模型，即將其簡化為X-Z平面上的二維等效模型。在此，僅施加ΔX和ΔZ兩個方向荷載的模型簡稱“2D”模型，而施加3個方向荷載的模型簡稱“3D”模型。

設(shè)管材屈服強(qiáng)度為235 MPa，管道壁厚為12 mm，外徑為1.824 m，埋深H為3.2 m，φ取14°，Su=66 kPa。僅改變地震震級進(jìn)行2D 和3D 模型的確定性分析結(jié)果對比。圖9（a）和（b）所示分別為不同地震震級作用下兩種模型所計算的管道最大拉應(yīng)變和最大壓應(yīng)變。由圖9可知，兩種模型計算所得管道最大拉應(yīng)變和最大壓應(yīng)變十分接近，說明了2D 模型的合理性。因此，在隨機(jī)模擬中采用2D模型替代3D模型是可行的。

3.3 結(jié)果分析

圖10（a）所示為根據(jù)1 000 組隨機(jī)樣本所得的管道最大拉伸應(yīng)變分布圖，其中橫軸表示管道最大拉伸應(yīng)變值，左側(cè)縱軸為每一個直方區(qū)間對應(yīng)的樣本頻率值，右側(cè)縱軸表示樣本累計頻率。采用2.4 節(jié)的應(yīng)變失效準(zhǔn)則判別，可得地震震級M 為6.3時，共有744 個失效樣本。圖10中隨機(jī)樣本的軸向最大拉應(yīng)變約有90%落在了0 到［εt］F之間，隨機(jī)樣本的軸向最大壓應(yīng)變有26.6%落在0 到［εc］F之間，即所有失效樣本均因局部受壓失效，且其中有10%的樣本既發(fā)生拉伸失效又發(fā)生壓縮失效。此外，在本文算例中，隨機(jī)變量均值的響應(yīng)（見2.5 節(jié)中的確定性分析）約為隨機(jī)模擬中響應(yīng)的中位數(shù)。當(dāng)發(fā)生6.3 級地震時，該管道的失效概率為74.4%，其主要是由于壓縮應(yīng)變較大引起，即發(fā)生壓縮破壞。這與確定性分析時的結(jié)果基本一致。在這種工況下，該管道不安全。

3.4 敏感性分析

除了考慮管道響應(yīng)分析中的各種不確定因素之外，本節(jié)探討地震震級、管壁厚度以及管道埋深對管道動力響應(yīng)的影響。如圖11所示，當(dāng)管壁厚度為12 mm，埋深H為3.2 m 時，隨著地震震級M的增大，管道的失效概率也急劇增大，管道的動力響應(yīng)對于地震震級十分敏感。

本文研究了地震震級M為6.3，埋深H為3.2 m，管壁厚度從10 mm 到20 mm 對管道動力響應(yīng)的影響。如圖12所示，隨著管壁厚度的增加，管道失效概率呈下降趨勢，在所選管道壁厚范圍內(nèi)管道的可靠性對于管道壁厚較敏感。

此外，本文研究了地震震級M為6.3 時，管道填土深度從1 m 到5 m 變化對管道動力響應(yīng)的影響。如圖13所示，隨著管道埋深的增加，未發(fā)現(xiàn)顯著的變化趨勢，與相關(guān)文獻(xiàn)［26-28］結(jié)論基本一致。由此可見，所研究的管道的安全性能對于管道埋深不敏感。

4 結(jié)果與討論

本文以滇中引水二期工程玉溪段管道五段為例，研究了跨斷層埋地管道的安全性和可靠性問題。采用ABAQUS 內(nèi)置管-土相互作用單元（PSI）建立了由地震位移控制的有限元模型，開展了地震作用下施加永久位移荷載時管道的響應(yīng)分析。考慮了地震位移荷載、管材強(qiáng)度和土體不排水剪切強(qiáng)度的不確定性，分析了管道動力響應(yīng)（軸向應(yīng)變）的分布，基于蒙特卡洛模擬計算了特定工況下管道的失效概率。此外，本文探討了地震震級、管壁厚度和填土深度對于所研究管道的失效概率的影響規(guī)律。主要結(jié)論如下。

（1）在本文算例中，忽略水平向位移ΔY（即將三維模型簡化為二維模型）對于跨斷層埋地管道響應(yīng)分析影響較小，主要是由于ΔY相對于ΔZ較小，且ΔY主要引起管道彎曲，對管道的最大軸向拉應(yīng)變和壓應(yīng)變的影響較小。因此，隨機(jī)模擬可采用2D模型，以提高計算效率。

（2）敏感性分析的結(jié)果表明，跨斷層埋地管道的動力響應(yīng)對于地震震級和管壁厚度較為敏感，管道的失效概率隨著地震震級的增加迅速增加。因此，在設(shè)計中應(yīng)合理確定地震震級。管道的失效概率隨著管壁厚度的增加逐漸下降，增加管壁厚度可以提高跨斷層管道的安全性和可靠性。另一方面，填土埋深對管道的失效概率的影響較小。

（3）在本文算例中，管道厚度增加到20 mm 時，管道的失效概率仍大于50%，管道的安全性和可靠性較低，僅依賴增加管壁厚度并不能有效提高管道的安全性和可靠性。本文算例中可以考慮采用強(qiáng)度更高的管材（比如Q345 或管線鋼等）或者減小斷層與管道的夾角以降低軸向位移荷載從而來提高工程的安全性。 □