基于離散元與有限元耦合的滑坡作用下輸氣管道動力響應分析

徐玲，彭星煜，魯廣航

(西南石油大學石油與天然氣工程學院，成都 610500)

油氣長輸管道具有輸送距離長、沿途地形地貌多樣、外部載荷復雜等特點[1-2]。在沿途地質災害破壞形式中，滑坡是主要形式也是破壞性極大的一種[3-4]?；聦τ蜌夤艿赖陌踩\行以及誘發(fā)管道失效造成的次生危害對自然環(huán)境和人類生命財產(chǎn)具有潛在威脅[5-6]，造成的經(jīng)濟損失不可控估量。因此能準確分析滑坡作用下油氣管道的動力響應對保障管道安全運行以及維修維護具有重要意義。

國內外學者針對滑坡災害對油氣管道的影響進行了一系列研究，目前研究手段主要包括解析解法、數(shù)值模擬以及物理模型三種。Philips等[7]、O’rourke等[8]、Chan[9]充分考慮了管道與滑坡土體間的相對滑動，運用土彈簧模型來表示管道與滑坡體的接觸作用，得到縱向、橫向兩種類型滑坡作用下管道受力分布情況。陳利瓊等[10]、王學龍[11]基于數(shù)值模擬分別采用CAESAR Ⅱ和ANSYS對埋地輸氣管道在縱穿和橫穿滑坡段下進行了應力分析與對比，并分析不同參數(shù)對管道應力應變的影響；朱勇等[12]、Daiyan等[13]通過物理實驗建立橫向滑坡模型，并結合彈性地基梁構建有限元模型，對比驗證了有限元模型的有效性，得出滑坡段管道應力分布情況。王榮有等[14]、王金安等[15]采用理論推導與數(shù)值模擬方法結合研究滑坡體的演化規(guī)律。當前大多數(shù)學者普遍采用數(shù)值模擬方法，管道與滑坡土體相互作用模型的選擇主要采用非線性接觸模型，但該模型未能充分考慮滑坡土體顆粒的離散性，不能準確反映出管道與滑坡體在接觸過程中的顆粒與管道、顆粒與顆粒之間的相互作用關系。模擬方法上采用單一的離散元或有限元方法，不能準確地模擬出顆粒與幾何結構間的接觸作用。

基于此，現(xiàn)提出采用離散元與有限元(discrete element method-finite element method，DEM-FEM)耦合分析法建立管道滑坡模型，從小尺度分析滑坡土體顆粒與管道的接觸作用，得到管道在滑坡作用下的受力情況，使模擬過程盡可能地符合工程實際，并通過具體工程案例進行研究分析。

1 離散元-有限元耦合

對于離散元與有限元的耦合方式，根據(jù)已有研究，主要有單向耦合和雙向耦合兩種。通過對比分析，選擇單項耦合方式，其示意圖如圖1所示。首先，利用離散元軟件模擬計算出顆粒與幾何結構之間產(chǎn)生的作用力；然后，將計算結果作為一種荷載條件導入到有限元軟件中；最后，結合其他約束條件對結構進行力學分析。這種方式忽略了幾何結構變形對顆粒的位置以及其與幾何體的接觸產(chǎn)生的影響，因此計算效率方面大大提高，聯(lián)合仿真也比較容易實現(xiàn)。而雙向耦合由于研究案例較少，在技術運用層面上還不夠成熟且操作過程較為復雜和困難，對計算機的性能要求極高，因此僅考慮單向耦合。

圖1 單向耦合示意圖

2 離散元基本方程

2.1 本構方程

離散元本構方程表征作用在兩實體上的接觸力與相對位移之間的關系?？偨佑|力由法向分量和切向分量組成。力與位移關系表達式為

(1)

2.2 運動方程

運動方程由平移運動方程和旋轉運動方程組成，表達式為

(2)

3 工程案例

聯(lián)絡線輸氣管道起于普光首站，經(jīng)平昌輸氣站，終于元壩站，管線全長202 km，管材為X70直縫埋弧焊鋼管，規(guī)格為Ф1 016 mm×17.5 mm，埋深為2 m，現(xiàn)實際運行壓力為7.5 MPa。由于省道409公路升級改造，在斜坡前緣進行開挖，加之施工期間當?shù)亟涤贻^多且較為集中，坡體產(chǎn)生變形滑移，造成管道滑坡事故，現(xiàn)場滑坡情況如圖2所示。

圖2 滑坡現(xiàn)場航拍圖

3.1 基本假設

在對實際工況進行模擬時，結合現(xiàn)場數(shù)據(jù)采集情況，需對模型做出合理的簡化與假設。

(1)在建立管道滑坡模型時，一般認為非滑坡區(qū)的寬度應大于滑坡區(qū)總寬度的1/4。

(2)假定滑移面為一平面，滑坡坡度為一定值，通過設置擋墻模擬邊坡，添加擋墻移動速度來模擬邊坡失穩(wěn)產(chǎn)生的滑坡效果。

(3)假定滑坡體內只含有土體顆粒而不含其他石塊等成分，滑坡過程僅考慮土顆粒對管道產(chǎn)生沖擊作用，顆粒形狀采用不規(guī)則的圓球組合模型。

(4)由于管道采用焊接方式連接，為降低計算難度，將焊縫的力學性能視為與母材一致。

3.2 滑坡計算模型建立

3.2.1 管道滑坡模型

根據(jù)滑坡輸氣管道現(xiàn)場情況，以及表1、表2基本模型參數(shù)，通過Solidwork建立滑坡輸氣管道三維模型，其俯視圖與側視圖如圖3所示。為保證管道在非滑坡區(qū)域內與基巖接觸充分，實現(xiàn)良好的嵌固，并綜合考慮仿真計算效率，建立管道模型總長為32 m，滑坡體寬度為10 m，基巖兩側的非滑坡區(qū)域分別為12 m。在滑坡體前人為設置擋墻模擬邊坡防護工程，通過設定擋墻移動速度來模擬因邊坡失穩(wěn)后顆粒因失去支撐而產(chǎn)生滑坡的工況。

表1 管道特性參數(shù)

表2 滑坡基巖與擋墻幾何參數(shù)

圖3 管道滑坡模型圖

3.2.2 土體顆粒模型

實際情況下滑坡土體顆粒形狀不規(guī)則，顆粒直徑大小也不盡相同，綜合考慮模型大小、計算效率等因素，采用簡單顆粒與復雜顆粒組合模型來模擬滑坡土體顆粒，顆粒模型如圖4所示，以便更加符合真實的滑坡體土體性質。顆粒參數(shù)參照文獻[16]的微觀參數(shù)標定實驗結果數(shù)據(jù)，如表3所示。采用系統(tǒng)中“Hertz-Mindlin with JKR”模型[17-18]來表征顆粒與幾何結構之間以及顆粒間的接觸特性。

圖4 顆粒模型圖

表3 顆粒特性參數(shù)

3.3 離散元模擬分析

3.3.1 離散元參數(shù)設置

由滑坡體幾何模型估算填充滿滑坡體所需顆粒，在EDEM軟件[19-20]中設置顆粒生成速率為20 000顆/s，生成的顆粒沿重力方向下落對滑坡體進行填充，顆粒填充過程總時長設為51 s。通過試算確定時間步長，在顆粒填充階段操作時間步設為Rayleigh時間步[21]的30%，在滑坡模擬階段操作時間步設為Rayleigh時間步的35%。數(shù)據(jù)存儲間隔在顆粒填充階段設置為1～2 s，在滑坡階段設置為0.5～1 s。通過虛擬顆粒工廠生成約102萬顆直徑為20 mm、位置均隨機分布的顆粒填充滑坡體。在此基礎上，設置擋墻移動速度為沿Z軸正方向3 m/s，滑坡過程模擬15 s，其他設置均保持不變，再選擇8線程提交運算。

3.3.2 離散元模擬結果

由離散元模擬得到滑坡過程中土體顆粒與管道不同時刻的接觸狀態(tài)如圖5所示，可以看出在滑坡發(fā)生3 s后管道前方土體顆粒與管道失去接觸，出現(xiàn)一定的臨空，此時管道后方土體無法及時移動至管道前方的空隙處，導致管道后部土體向下滑移并擠壓管道。隨著滑坡的繼續(xù)，管道前方的臨空面逐漸增大，并且逐漸向重力方向偏轉，而管道后部土體顆粒因管道外表面的阻擋無法繼續(xù)滑移，持續(xù)推擠管道，直到滑坡結束，管道前方外表面已經(jīng)完全臨空。

圖5 滑坡過程中不同時刻顆粒-管道三維接觸狀態(tài)

通過EDEM后處理得到管道在滑坡過程中受到的合力以及各方向受力時程變化曲線如圖6所示，可以看出在滑坡發(fā)生后管道所受合力出現(xiàn)先增大后減少的趨勢，管道在軸向(X方向)上受力較小，可忽略其影響，在重力方向(-Y方向)上由于管道上方的顆粒的減少，管道受到該方向上的作用力隨時間逐漸減小，直至管道部分暴露出埋土外，管道受力保持穩(wěn)定。此外，管道在迎坡方向(Z方向)上受力變化趨勢與管道所受合力變化趨勢一致。由圖可以看出，管道所受最大合力與迎坡方向受到的最大作用力均發(fā)生在第3秒，最大合力為 53 686.7 N，重力方向最大受力為53 686.7 N，迎坡方向最大受力為33 822.7 N，管道在滑坡作用下所受最大推力方向與重力方向夾角為36°。

圖6 滑坡過程中管道受力時程曲線

3.4 有限元模擬分析

3.4.1 本構模型與參數(shù)設置

基巖本構模型選擇目前普遍使用的Mohr-Coulomb模型，管材本構模型采用Ramberg-Osgood模型。設置基巖與管道間為摩擦接觸，根據(jù)研究摩擦因數(shù)設為0.3，將管道在第3秒時刻各單元受到的最大滑坡推力結果文件通過插件導入ANSYS，軟件自動識別出滑坡推力作用在管道上的空間分布，管道內壓設置為7.5 MPa，并且考慮重力場的作用，對基巖底部添加固定約束，其作用面選擇滑坡區(qū)域管道外表面，完成求解設置后提交運算。

3.4.2 有限元模擬結果

在滑坡作用下，由于滑坡土體顆粒對管道的沖擊，由圖7(a)可以看出，在滑坡體主滑線中部受到最大等效應力為237.53 MPa；在滑坡過程中，滑坡體與基巖在軸向會形成相對運動，在滑坡邊界基巖會對管道產(chǎn)生較大的剪切作用，因此在該處管道受到剪切應力最大，由圖7(b)可得最大剪切應力為60.555 MPa；同時由圖7(c)可以看出，管道在滑坡體主滑線附近的變形位移最大為12.495 mm，處于非滑坡區(qū)域內管道變形位移為零，由此表明非滑坡區(qū)管道在基巖中處于嵌固狀態(tài)，因此可以說明兩側非滑坡區(qū)管道長度設為11 m長能夠有效實現(xiàn)嵌固作用。由圖7(d)滑坡段管道的等效應力分布可知，在滑坡管段中部和滑坡周界處的應力最大。

圖7 管道應力分析

3.5 模型驗證

通過離散元與有限元的耦合得到滑坡過程中管道受到的最大滑坡推力與應力分布情況，結合管道應力監(jiān)測平臺，平臺設有3個監(jiān)測點，分別位于滑坡兩邊界處(監(jiān)測點1、3)和滑坡管道中部(監(jiān)測點2)。對比3處管道2020年4月4—9日的應力監(jiān)測數(shù)據(jù)，如圖8所示，在同一滑坡推力下，滑坡管道中部受到的應力最大監(jiān)測值為224.52 MPa，滑坡邊界處的應力最大監(jiān)測值分別為277.91 MPa和286.84 MPa；有限元模擬值分別為237.53 MPa和299.79 MPa(左右邊界對稱)。由此可見，現(xiàn)場監(jiān)測數(shù)據(jù)和有限元模擬結果得到的應力大小趨勢基本一致，且3個監(jiān)測點的相對誤差在7.30%以內，在可接受的范圍中，因此驗證該模型具有一定的可靠性和參考意義。

圖8 監(jiān)測平臺實時監(jiān)測數(shù)據(jù)

4 結論

綜上分析，可得出以下結論。

(1)通過離散元和有限元單向耦合方式，基于EDEM模擬土體顆粒滑坡過程，得到管道在各方向的受力時程變化曲線，將第3秒時刻管道所受滑坡最大推力作為載荷導入ANSYS中，模擬得到管道等效應力、剪切應力以及位移的分布云圖。應力在滑坡管段中部和滑坡周界處的最大，變形位移在滑坡體主滑線附近最大。

(2)將數(shù)值模擬得到的滑坡管道中部和滑坡周界處的最大等效應力值與對應位置的監(jiān)測點數(shù)據(jù)進行對比，相對誤差在7.30%以內，因此表明模型具有一定的可靠性。

(3)根據(jù)滑坡段管道等效應力分布，輸氣管道若不可繞避滑坡體，則建議加強對滑坡管道中部和滑坡周界處，以及環(huán)焊縫連接處的實時監(jiān)測，避免應力集中造成管道加速破壞。