基于非線性有限元和響應面方法的海底管道受墜物撞擊損傷的風險分析

姜逢源， 董 勝， 趙玉良

(中國海洋大學 工程學院，山東 青島 266100)

海底管道被視為近海油氣產(chǎn)業(yè)的“生命線”，其作為連接海底設施與陸上終端的介質(zhì)，承擔著輸送油氣的重要任務，具有快捷、高效等優(yōu)點。由于人類海洋活動越發(fā)頻繁，平臺作業(yè)和船舶拋錨等第三方活動容易對海域內(nèi)的管道造成墜擊損傷，從而產(chǎn)生油氣泄露，這已成為引發(fā)管道失效事故的主要原因[1]。由于海底管道所處環(huán)境復雜，檢修困難，一旦失效，將帶來不可估量的經(jīng)濟損失和環(huán)境破壞。因此，如何確定管道正常營運狀態(tài)下的可靠度，準確評估管道的失效風險是一項重要工作。

迄今，針對第三方破壞引發(fā)的海底管道失效風險已有一定的研究。有學者對于墜物落水及船舶撞擊管道的概率給出了確定方法[2-3]。Bai等[4]估計了多種類型撞擊下海底管道的失效概率并研究其風險預測及可接受準則。Mazzola[5]提出了一種計算海底管道受平臺墜物撞擊的失效概率的方法，分析了不同墜落點、管道地點對失效概率的影響。DNV[6]提出了墜物撞擊下海底管道失效概率的計算方法(以下稱DNV法)，該方法的提出為后續(xù)研究提供了計算框架，由于采用經(jīng)驗公式計算管道的損傷，忽略了諸多非線性因素的影響。Kawsar等[7]采用隨機抽樣法模擬了不同工況，結(jié)合有限元法求得管道失效的超越概率，但未考慮管土相互作用的影響。丁紅巖等[8]和Jiang等[9]對DNV法進行了改進，引入可靠度理論求解失效概率，分析了不同墜物對管道失效概率的影響，但二者對管道失效狀態(tài)的模擬均采用計算公式法，忽略了非線性因素的影響。

受海底環(huán)境荷載作用和邊界條件的影響，海底管道在沖擊荷載下的響應涉及到的非線性作用復雜，需通過非線性有限元分析確定其失效狀態(tài)[10-11]。目前，通常采用經(jīng)驗公式法來描述海底管道的結(jié)構(gòu)失效狀態(tài)[12-13]，由于忽略了諸多非線性因素，因而無法反映結(jié)構(gòu)的真實狀態(tài)。海床土體作為支撐管道的介質(zhì)，對沖擊能力的耗散起到可觀的作用[14]，目前相關分析中，對上述管-土耦合作用的考慮相對欠缺，計算得到的管道損傷較實際情況偏大，導致分析結(jié)果較為保守，造成不必要的資金投入。其中，挖溝埋深作為保障管道安全的重要手段，國內(nèi)外行業(yè)規(guī)范中對于埋深的規(guī)定[15]始終沒有明確標準，至今尚未建立起與失效概率與埋深的定量關系，無法有效地指導工程設計。

綜上所述，目前關于管道損傷的風險分析方法對非線性作用考慮不充分，對管土耦合作用的考慮更是寥寥無幾，影響工程設計與風險控制。為解決以上問題，基于Python語言的二次開發(fā)技術[16]，將非線性有限元分析法與響應面法[17]耦合，建立海底管道受墜物撞擊損傷的可靠度分析模型，探討了非線性因素及管-土耦合作用對管道失效概率的影響。本文首先介紹了可靠度分析模型的建立，其次將有限元計算結(jié)果及失效概率估計值分別與物理模型試驗值及Monte Carlo法的估計值對比，驗證了該模型合理性。隨后從工程應用的角度，考慮管-土耦合作用因素，分析不同條件下的失效概率。最后對隨機變量進行敏感性分析，得出結(jié)論。

1 非線性有限元分析

1.1 顯示動力學算法

海底管道的撞擊屬于高速動力學問題，荷載持續(xù)時間短，量級大，結(jié)構(gòu)反應變化快，涉及高度非線性因素。顯式動力學算法可精確捕捉極小時間步內(nèi)結(jié)構(gòu)的動力響應，適用于此類問題的求解。結(jié)構(gòu)的控制方程如式(1)所示，將其中的加速度項及速度項采用中心差分代替(見式(2)和式(3))，可得到相應的線性方程組(見式(4))，對其求解得到單元節(jié)點位移U后，將其代入幾何及物理方程中，即可求得應力及應變。

(1)

(2)

(3)

(4)

其中，

(5)

(6)

式中：M，C，K，Rt分別為系統(tǒng)的質(zhì)量矩陣、阻尼矩陣、剛度矩陣、荷載矩陣；Ut為單元節(jié)點位移矩陣。

1.2 耦合歐拉-拉格朗日算法

墜物侵入管道上覆土體過程中涉及高度非線性、多介質(zhì)接觸等問題，使土體產(chǎn)生極大變形，傳統(tǒng)有限元分析方法會產(chǎn)生網(wǎng)格畸變，難以收斂。耦合歐拉-拉格朗日算法[18](coupled Eulerian-Lagrangian, CEL)是目前結(jié)構(gòu)物-土體相互作用分析中的常用方法，且已嵌入有非線性限元分析程序Abaqus中。其利用歐拉網(wǎng)格在空間保持固定，土體材料在其中自由流動的建模方式有效解決土體大變形問題，并以歐拉體積分數(shù)EVF追蹤材料的位置。同時采用拉格朗日網(wǎng)格去描述結(jié)構(gòu)物，其應力應變響應則通過歐拉-拉格朗日接觸算法獲取。如圖1所示，通過計算歐拉材料節(jié)點和拉格朗日節(jié)點在每一時間步內(nèi)的相對位移來確定作用在二者之間的罰力，完成相互間信息的映射傳遞。其計算公式為

Fp=kp·dp

(7)

(8)

FE,i=NiβiFp

(9)

(10)

圖1 耦合歐拉-拉格朗日接觸算法原理Fig.1 Contact algorithm principle of CEL method

基于Abaqus平臺，本文建立的有限元模型如圖2所示，管道采用拉格朗日網(wǎng)格進行離散，其鋼材料的應力應變關系使用雙線性各向同性強化模型進行描述；海床土體則用歐拉網(wǎng)格離散，并對海床土體施加初始地應力及海水壓力場，土體的彈塑形行為則用Mohr-Coulomb(MC)本構(gòu)模型表示；對墜物同樣采用拉格朗日網(wǎng)格離散，施加剛體約束及初速度，并以等效密度的方式考慮墜物的水動力因素。

注： D為管徑； e為埋深。圖2 墜物-管道-土體相互作用有限元模型Fig.2 Finite element model of interactions between dropped object-pipeline-soil

1.3 有限元模型驗證

為驗證非線性有限元模型的可靠性，選取文獻[19]中船錨與埋深管道撞擊的物理模型試驗進行對比驗證，其結(jié)果如圖3所示。其中試驗土樣為黏土，埋深e由0～0.2 m變化。由對比結(jié)果可知，總體上數(shù)值模型與試驗的結(jié)果吻合較好，平均相對誤差為6.58%。但尚存在些許誤差，如埋深e=0.1 m時的情況。其原因很可能是由于不同埋深管道周圍土體的固結(jié)狀態(tài)難以保持一致，使得土體力學特性的真實值與測量值之間存在差異，從而使得計算結(jié)果與試驗結(jié)果產(chǎn)生偏差。由上述分析可知，本文的有限元模型可較為精確的對管道損傷進行預測。

圖3 管道凹痕深度對比Fig.3 The comparison of pipeline dent depths

2 有限元-響應面耦合可靠度模型

2.1 功能函數(shù)

DNV規(guī)范中，海底管道的凹痕深度δ與管徑D的比值是評估結(jié)構(gòu)是否失效的重要指標，因此管道的極限狀態(tài)方程可表示為式(11)，其中η為無量綱臨界參數(shù)，用以控制失效判定條件，在此取η=0.2以考慮管道發(fā)生破裂造成較大危害的情況。

Z=g(x)=δ(x)/D-η

(11)

式中：x為隨機變量；凹痕深度δ(·)為隨機變量的函數(shù)，經(jīng)有限元分析獲得。

由于海底條件復雜，海底管道的制造、加工及施工過程較為復雜，難免出現(xiàn)一些偏差，從而影響管道的承載能力及幾何尺寸，使上述因素存在變異性。故在本文的可靠度分析中，將管道鋼材料的屈服強度σy、管徑D及壁厚t作為隨機變量來處理。上述三個隨機變量均服從正態(tài)分布，均值取其標準值，變異系數(shù)依次為0.07，0.03，0.05。

2.2 實現(xiàn)步驟

由以上敘述可知，海底管道在沖擊荷載下的損傷涉及復雜的非線性因素，結(jié)構(gòu)的極限狀態(tài)方程無法顯式表達，傳統(tǒng)的一次二階矩法難以進行求解。針對極限狀態(tài)方程為隱式的情況，常用的方法有：Monte-Carlo法[20]、神經(jīng)網(wǎng)絡法[21]和響應面法。但前兩種方法均需大量的抽樣計算，不適用于大型計算模型的求解。而響應面法僅需在結(jié)構(gòu)的失效域附近進行一系列數(shù)值計算，構(gòu)造出一個具有顯式表達式的函數(shù)來近似結(jié)代替結(jié)構(gòu)真實的極限狀態(tài)曲面，從而可利用一次二階矩法對可靠度進行求解，具有明確的收斂準則和較高的計算效率。

(12)

(13)

Python語言作為Abaqus二次開發(fā)支持的腳本語言，可用于有限元分析與響應面法的耦合，完成有限元模型的建立與更新、計算結(jié)果的提取、求解可靠度設計驗算點、中心點的搜索及循環(huán)計算，使可靠度求解程序化，無需進行人為干預。表1給出了Python語言的實現(xiàn)流程。

表1 可靠度分析模型的Python實現(xiàn)

2.3 可靠度模型驗證

為驗證本文提出的可靠度分析模型的合理性，選取表2所列出的12種工況進行可靠度指標的計算，其中為便于對比，海床土體簡化為剛性海床，墜物簡化為球狀剛體。將計算結(jié)果與Monte Carlo法進行對比。由于Monte Carlo法需進行大量抽樣計算，考慮到時間成本及可行性，該方法中對于管道凹痕采用DNV規(guī)范推薦的公式[23]進行計算

(14)

式中：E為墜物的沖擊能量；D為管道直徑；t為管道壁厚；mp=1/4(σy·t2)為塑性彎矩，σy為管道鋼材料屈服強度。

表2 工況設計

計算結(jié)果如圖4所示，大多數(shù)計算點分布在y=x直線上部附近，這主要是由于DNV推薦的理論公式假設墜物為尖銳的楔體，因而計算得到的凹痕值較大，導致Monte Carlo法計算的可靠度指標偏小。絕大多數(shù)計算點落在±20%的相對誤差線范圍內(nèi)，僅有2個點落在誤差線外，但其可靠度指標的絕對誤差分別為0.05和0.06。所有工況中絕對誤差最大為0.500，最小為0.05，平均為0.217；相對誤差最大為21.39%，最小為4.23%，平均為13.67%。

圖4 可靠度指標對比Fig.4 The comparisons of reliability index

分析誤差產(chǎn)生的原因，可歸納為兩類因素：①來自于管道凹痕深度計算方法的差異。本文提出的方法中采用非線性有限元法計算管道的損傷，可考慮材料、接觸條件及邊界條件的非線性因素，因而可獲得更為精確的結(jié)果。而Monte Carlo法采用的則是計算公式法，對實際條件進行了簡化，忽略了墜物與管道間的相互作用，無法客觀描述物體間的接觸，尤其當需考慮墜物形狀及管道-海床相互作用時會與實際情況存在較大偏差。②可靠度算法的差異對計算結(jié)果造成一些影響。響應面法主要在管道的失效域附近設置試驗點，有些情況可能會影響設計驗算點的收斂效果從而引起部分偏差。

從總體上看，該可靠度計算模型的結(jié)果與Monte-Carlo法的結(jié)果相差不大，平均相對誤差控制在15%以內(nèi)，說明該模型的計算結(jié)果是合理的，同時考慮了材料、接觸等非線性因素，更為符合實際情況。

3 分析討論

3.1 管土耦合作用對損傷的影響

海底管道通常安置在海床表面或埋置于海床土體中，現(xiàn)分析上述因素對管道損傷的影響?？紤]海床土體為均質(zhì)黏土，計算管道在柔性海床及埋深(e=0.5 m)條件下受墜物撞擊時(E=28.9 kJ)的損傷情況，并與剛性海床條件下的情況進行對比。

計算結(jié)果如圖5所示，從兩方面進行分析。①管道損傷的分布：在剛性海床的條件下，由管道的位移云圖可看出，管道損傷的凹痕深度在撞擊中心處最大，并向四周擴散，近似成橢圓形分布，此時沖擊能量幾乎全部被管道塑性變形所吸收。而在海床柔性及埋深的條件下，位移云圖的分布范圍逐漸向管道撞擊中心靠攏，其損傷區(qū)域的面積不斷減小。②管道損傷的最大凹痕深度：管道裸落在剛性海床上變?yōu)槁裰糜诤４矁?nèi)，撞擊中心處橫截面的凹痕深度由則0.18D逐漸減小至0.06D。這是由于在撞擊過程中，一部分沖擊能量被管道周圍的土體變形及管道的整體位移所耗散，減緩了管道的損傷范圍及深度，這與Zeinoddini等和姜逢源等得到的結(jié)論是一致的。分析可知，海床柔性及埋深對管道的損傷影響較大，不考慮該因素時得到的結(jié)果偏保守。因此，在可靠度分析中需考慮上述因素：使用Python編程建立相應條件下的有限元模型并嵌入到可靠度分析程序。

圖5 不同條件下管道的位移云圖及凹痕深度Fig.5 Pipeline displacement contour and dent depth in different conditions

3.2 實例分析

圖6為某管道工程的平面布置圖，該工程所在地水深d=100 m，所在海域墜物的平均質(zhì)量為900 kg左右；海床土體為均質(zhì)黏土(不排水抗剪強度Su=10 kPa)；管道西南側(cè)距墜物點的距離為60 m，向東北方向延伸150 m，管道鋼材料的屈服強度σy=360 MPa,管道外徑D=324 mm，壁厚t=12 mm。

圖6 海底管道布置圖(m)Fig.6 Filed layout of submarine pipeline(m)

現(xiàn)分析管土相互作用對管道失效風險的影響，簡化墜物為質(zhì)量m=900 kg、速度v=10 m/s球狀墜物，計算該管道在三種海床條件(①剛性海床；②柔性海床；③埋深e=0.5 m)下，受墜物撞擊時的可靠度指標β及相應的失效概率Pf=1-Φ(β)。結(jié)果如圖7所示。三種條件下管道的可靠度指標分別為-4.85，2.93和5.21。海床柔性及埋深會大大提高管道的可靠度指標，使其失效概率由1逐漸降低至10-7，對海底管道進行埋深處理是一種有效的保護措施。

圖7 不同海床條件下管道的可靠度指標及失效概率Fig.7 Reliability index and failure probability for pipeline under different seabed conditions

根據(jù)圖6中管道的鋪設路徑及墜物點的空間位置關系可求出墜物撞擊管道的概率Ph，具體計算方法可參考DNV法，則最終管道的總失效概率為Pf,t=Ph×Pf。圖8給出了三種條件下管道總失效概率沿管道長度方向的空間分布，并與DNV法的結(jié)果進行了對比。由結(jié)果可知：①剛性海床條件下管道總失效概率的計算結(jié)果與DNV法的結(jié)果十分接近，因為二者均忽略了管土相互作用的影響；②失效概率隨著管道距起點距離的增加迅速減小，對于剛性海床情況，距起點的距離大于40 m后其失效概率小于10-5，滿足一般工程的安全標準，故應對前40 m內(nèi)的管道段進行重點保護；③柔性海床條件下最大的失效概率為10-6，為剛性海床條件下的千分之一，可見DNV法未考慮土體對沖擊能量的耗散，計算結(jié)果偏于保守；④埋深條件下，管道的最大失效概率僅為10-10，這是因為管道上覆土體的剪切破壞過程需消耗大量能量，最終有效保護了管道。

3.3 埋深對失效概率的影響

由3.2節(jié)分析可知，埋深極大提高了管道可靠性。為進一步探究其影響，另選3種工況，計算不同埋深條件下(e=0.5 m，0.75 m，1.0 m)管道的可靠度指標及失效概率，計算結(jié)果如圖9所示。其中，墜物為球狀，相關參數(shù)為：①m=3 500 kg，v=9.06 m/s；②m=2 750 kg，v=8.71 m/s；③m=2 000 kg，v=8.07 m/s。管道參數(shù)：σy=538 MPa，D=324 mm，t=12 mm?？梢钥闯?，管道受撞擊時的失效概率隨著埋深的增加而呈非線性降低。對于工況2的情況，0.5 m或0.75 m的埋深時管道的失效概率均大于10-5，其不足以有效的保護管道，直至埋深增加至1 m時，其失效概率才控制在可接受的范圍內(nèi)。但對于工況1來說，沖擊能量進一步增大，1 m的埋深也不足以保護管道，還需增加埋深或采取其他的保護措施。因此，為選擇合理的管道埋深，應以該海域內(nèi)典型的墜物沖擊為目標進行設計。

圖8 管道失效概率沿管道長度的空間分布Fig.8 Spatial distribution of failure probability along the length of pipeline

圖9 不同埋深下管道的可靠度指標及失效概率Fig.9 Reliability index and failure probability for pipeline with different burial depth

此外，Zeinoddini等和姜逢源等研究中還指出海床土體的力學特性也對沖擊能量的耗散有著重要的影響，故對其也應進行考慮。圖10為埋深e=0.5 m時，在工況2條件下管道的可靠度指標及失效概率隨土體強度的變化曲線?？梢钥闯?，隨著土體強度的增加，管道的失效概率逐步下降，并且下降的趨勢隨著土體強度的增加而越發(fā)明顯，這可從姜逢源等研究中的結(jié)果得到印證：對于埋置管道，隨著土體強度的增加，管道損傷的凹痕值的減少量逐漸增大。

3.4 隨機變量的靈敏度分析

進行可靠度分析時，選取了σy，D，t作為隨機變量。由于設計和施工等因素，相關參數(shù)需要進行調(diào)整，為分析上述因素的變異性的影響，以3.3節(jié)中工況3為例，即以μσy=448 MPa，cvσy=0.07；μD=324 mm，cvD=0.03;μt=12 mm，cvt=0.05為基準，分別取其0.8倍、0.9倍、1.0倍、1.1倍、1.2倍，計算在質(zhì)量為m=2 000 kg的球狀墜物撞擊(v=8.07 m/s)埋深e=0.5 m的海底管道時的可靠度及失效概率，結(jié)果如圖11和圖12所示。

圖10 管道的可靠度指標及失效概率隨土體強度的變化Fig.10 Reliability index and failure probability for pipeline with different soil strength

圖11 可靠度指標及失效概率隨均值的變化Fig.11 Reliability index and failure probability with different means

分析可得：①當均值發(fā)生變化時，屈服強度σy的影響效果最為顯著，可靠度指標從1.85上升至4.28，其次為壁厚t，再次為管徑D。說明管道材料的屈服強度的增加明顯提高了管道抵抗沖擊的承載力，而壁厚與塑性彎矩呈2次方關系，其影響效果也較為明顯。②當變異系數(shù)發(fā)生變化時，可靠度指標隨著cvσy的增加而下降，其失效概率發(fā)生了量級上的變化，由10-4迅速增加到10-2，因而在對鋼材的力學強度進行檢測時要嚴格把控其質(zhì)量。而對于cvD及cvt的變化，可靠度指標僅發(fā)生微小的波動，受管徑D及壁厚t變異性的影響較小。

圖12 可靠度指標及失效概率隨均值的變化Fig. 12 Reliability index and failure probability with different variable coefficients

4 結(jié) 論

本文將非線性及管土耦合作用引入風險分析中，利用非線性有限元分析法與響應面法耦合，構(gòu)建了海底管道受墜物撞擊的可靠度分析方法，并與Monte Carlo法對比，驗證了其合理性，同時分析了管土耦合作用對管道失效概率的影響，得到以下結(jié)論：

(1)基于Python二次開發(fā)技術，將非線性有限元分析模塊嵌入響應面法中構(gòu)建了可靠度分析模型，實現(xiàn)了求解算法與有限元分析之間的信息傳遞、迭代等工作，該技術可拓展到其他相關領域，如結(jié)構(gòu)的參數(shù)優(yōu)化等。

(2) 由實例分析可知，管道的失效概率沿管道起點距離逐漸減小，前40 m范圍內(nèi)的失效概率較大，應重點進行保護；與剛性海床條件相比，柔性海床條件下管道的失效概率為前者的10-3；埋深條件下，管道的最大失效概率僅為剛性海床條件下的10-7，極大提高管道的安全性；綜上，不考慮管土相互作用時結(jié)果偏保守；

(3) 增加埋深能明顯降低管道的失效風險，但沖擊能量較大時(約145 kJ)，1 m埋深時管道的失效概率仍大于10-5，不足以提供安全保障。故在選擇合理的埋深時，應針對該海域內(nèi)墜物的典型墜物附帶的沖擊能量作為設計埋深，同時也要考慮到海床土體的力學特性。

(4) 隨機變量敏感性分析，當各變量的均值以0.8～1.2倍變化時，σy對可靠度及失效概率的影響最大，其次為t及D；當變異系數(shù)以0.8～1.2倍變化時，可靠度指標隨cvσy增大而減小，失效概率發(fā)生量級上的增加；其余兩變量的變異系數(shù)對結(jié)果的影響不大，工程建設中，檢測時應嚴格把控鋼材的材質(zhì)。

本文提出的有限元-響應面耦合可靠度模型首次應用于管道的風險分析中，建立了管道風險與埋深間的關系，該模型及研究結(jié)果可為工程設計提供有效指導。未來，將在模型中考慮更多因素的作用，如墜物種類、尺寸，土體材料等影響。