基于離散元與有限元耦合的滑坡作用下輸氣管道動(dòng)力響應(yīng)分析

徐玲，彭星煜，魯廣航

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

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

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

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

1 離散元-有限元耦合

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

圖1 單向耦合示意圖

2 離散元基本方程

2.1 本構(gòu)方程

離散元本構(gòu)方程表征作用在兩實(shí)體上的接觸力與相對(duì)位移之間的關(guān)系?？偨佑|力由法向分量和切向分量組成。力與位移關(guān)系表達(dá)式為

(1)

2.2 運(yùn)動(dòng)方程

運(yùn)動(dòng)方程由平移運(yùn)動(dòng)方程和旋轉(zhuǎn)運(yùn)動(dòng)方程組成，表達(dá)式為

(2)

3 工程案例

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

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

3.1 基本假設(shè)

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

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

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

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

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

3.2 滑坡計(jì)算模型建立

3.2.1 管道滑坡模型

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

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

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

圖3 管道滑坡模型圖

3.2.2 土體顆粒模型

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

圖4 顆粒模型圖

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

3.3 離散元模擬分析

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

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

3.3.2 離散元模擬結(jié)果

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

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

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

圖6 滑坡過(guò)程中管道受力時(shí)程曲線(xiàn)

3.4 有限元模擬分析

3.4.1 本構(gòu)模型與參數(shù)設(shè)置

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

3.4.2 有限元模擬結(jié)果

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

圖7 管道應(yīng)力分析

3.5 模型驗(yàn)證

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

圖8 監(jiān)測(cè)平臺(tái)實(shí)時(shí)監(jiān)測(cè)數(shù)據(jù)

4 結(jié)論

綜上分析，可得出以下結(jié)論。

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

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

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