埋地管道三參量粘彈性邊界地震模型及其動力響應研究

吳靜波，杜 強，曾祥國，李 偉，趙清龍

(1.四川大學 建筑與環境學院，深地科學與工程教育部重點實驗室， 成都 610065；2.西南油氣田分公司通信與信息技術中心， 成都 610051)

埋地管道是天然氣運輸的主要途徑，具有運量大、占地少的優點。埋地管道分布廣泛，其承受的外界載荷錯綜復雜，國內外天然氣管道事故時而發生[1-3]，地震也是造成輸氣管道破壞的重要因素之一[4]。地震載荷作用下，埋地輸氣管道一旦被破壞，將嚴重影響能源的運輸和供應，同時逸出的天然氣還可能引起火災和爆炸，造成巨大經濟損失。因此，研究埋地輸氣管道設施在地震荷載作用下的動力響應并進行安全性評價十分必要。

地震作用時，結構瞬間承受過大的載荷，極易引發結構的失效和破壞。由于地表結構的破壞和損傷比地下結構更容易觀測到，現有的研究更多集中在地震對地表結構的影響，如地面建筑物、大壩等。然而，在近些年的一些地震[5-7]中，出現了嚴重的地下結構破壞現象，這促使地震工程師們更多關注地下結構的設計[8-12]。在進行地震動力分析時，結構與無限域的地基存在能量交換，模型邊界的處理對計算結果影響很大。Deeks等[13]和劉晶波[14-15]建立了粘彈性人工邊界，在模擬結構-地基的動力響應時具有較高的精度。

本文在原有的粘彈性邊界元件上串聯一個彈簧元件，建立了三參量粘彈性人工邊界，使用商業有限元軟件ANSYS實現并驗證了三參量粘彈性邊界。最后，將三參量粘彈性邊界應用到埋地管道的地震模型中，研究了地震波入射角和管道埋深對埋地管道動力響應的影響，為埋地輸氣管道的抗震設計提供了一定依據。

1 三參量粘彈性人工邊界

1.1 現有的粘彈性人工邊界及其缺點

如圖1所示，現有地震模型使用的粘彈性人工邊界普遍采用兩參數Kelvin固體模型。在進行動力分析時，將管道和周邊一定區域內的土壤從無限域中切取出來，并在邊界的法向(圖1中的Y向)和2個切向(X向和Z向)施加彈簧阻尼器系統來連接有限的管道結構和地基，使邊界上滿足相當于無限域的應力-位移協調條件。邊界上彈簧阻尼器系統的彈簧剛度和阻尼系數的取值為[16-17]：

(1)

(2)

圖1 管道結構的兩參數粘彈性邊界地震模型示意圖

上述的兩參數粘彈性人工邊界考慮了對散射波能量的吸收和邊界外半無限介質的彈性恢復，在土-結構相互作用模型的地震分析中得到了廣泛應用。然而，這種模型存在致命的弱點，兩參數粘彈性邊界采用的Kelvin固體模型無法描述松弛特性[19]，即不能描述巖土完整的蠕變、松弛、回復的粘彈性特性。

1.2 三參量粘彈性人工邊界地震模型的建立

在同時考慮人工邊界對散射波能量吸收和邊界外半無限介質彈性恢復的前提下，為了能完整描述邊界層的粘彈性，本文在原有的粘彈性邊界元件上串聯一個元件，構成了三參量固體(Kelvin-Voigt模型)模型(該模型能描述完整的蠕變、松弛、回復的粘彈性特性[19])，并使用該模型建立三參量粘彈性人工邊界。

1.2.1有限地震模型法向邊界等效參數的推導

將結構從無限域介質中截斷，并在底部和四周的邊界節點上施加三參量粘彈性邊界。在邊界的法向，施加如圖2所示的彈簧-彈簧-阻尼器-集中質量系統。該物理系統的運動平衡方程和位移協調方程為：

(3)

(4)

式中:U為邊界節點的位移；σ為位移U在人工邊界處產生的應力；K1和K2分別為圖2中對應位置彈簧的剛度；U1為與邊界節點相鄰彈簧的位移；U2為遠離邊界節點彈簧的位移；UC為阻尼器的位移；UM為集中質量點的位移；C為阻尼器的阻尼系數；M為集中質量。

圖2 邊界法向的三參數粘彈性物理系統示意圖

由式(3)和(4)可得邊界節點應力和位移的關系：

(5)

(6)

球坐標系中，P波波陣面上法向應力與位移的關系[16]為：

(7)

式中：R為人工邊界節點與散射源的距離；G為介質的剪切模量；ρ為介質的密度；cP為介質的P波波速。

對照式(6)和(7)，有：

(8)

1.2.2有限地震模型切向邊界等效參數的推導

對于邊界的切向，直接施施加三參量固體(Kelvin-Voigt模型)模型即可，如圖3。

圖3 邊界切向的三參數粘彈性物理系統

Kelvin-Voigt模型應力與位移的關系為：

(9)

(10)

球坐標系中，S波波陣面上切向應力與位移的關系[16]為：

(11)

式中：cS為介質的S波波速。

對比式(10)和(11)，可得：

(12)

1.2.3數值模擬技術中的參數設置

圖2中的質量M和阻尼器組成一個不穩定的系統，為了方便計算和處理，數值模擬中將質量M忽略，并將阻尼器固定。三參量粘彈性邊界可以通過有限元軟件ANSYS的Combin14單元實現。具體實現方法為：在需要施加人工邊界的截面上向法向和切向均延伸兩層Combin14單元，其中一層單元設置同時設置彈簧剛度和阻尼系數，另一層單元只設置彈簧剛度。法向的彈簧剛度KN1、KN2、阻尼系數CN和切向的彈簧剛度KT1、KT2、阻尼系數CT分別為：

(13)

(14)

1.2.4三參數粘彈性邊界的地震輸入方式

數值模擬時，將地震波加速度轉換為人工粘彈性邊界上的等效節點載荷來實現地震波的輸入。在邊界上應力和位移滿足圖4所示的關系：

(15)

(16)

式中：FB為施加在人工邊界上的等效節點載荷；fB為彈簧阻尼器系統與人工邊界連接處的內力；U為地震波在人工邊界處產生的位移；σ為位移U在人工邊界處產生的應力；U1為與邊界點相鄰彈簧的位移；U2為遠離邊界點彈簧的位移，可以由位移U求得。

于是，人工粘彈性邊界上的等效節點載荷FB為：

(17)

圖4 三參量人工粘彈性邊界上的等效載荷

1.3三參量固體(Kelvin-Voigt模型)粘彈性人工邊界的驗證

本文使用ANSYS的combin14單元，并用ANSYS命令流語言編寫程序實現了三參量粘彈性人工邊界及其對應的地震波輸入。然后，使用有限元軟件ANSYS模擬了三維Lamb表面源問題和P波斜入射的波動問題，得到的位移數值解與(近似)理論解吻合得較好，驗證了三參量粘彈性邊界的有效性。

1.3.1三維Lamb表面源問題

三維表面源問題的數值模型如圖5所示，模型范圍為X向-60～60 m，Y向-60～0 m，Z向-60～60 m。介質為均勻線彈性材料，其剪切模量為500 MPa，泊松比為0.25，密度為2 000 kg/m3。模型的上表面為自由面，底部和4個側面考慮了三參量粘彈性邊界(如圖5，在模型邊界節點的法向和2個切向共施加3組3個元件的彈簧阻尼器系統)、兩參量粘彈性邊界與遠置邊界(作為近似解)。 三參量粘彈性邊界的n取10 000，法向、切向邊界參數取4/3和2/3；兩參量邊界的法向、切向邊界參數取4/3和2/3；遠置邊界將原有邊界向外延伸了500 m(在計算時間 1 s內，反射波未傳回計算區域)。在坐標原點施加動力載荷F(t)，并觀測點A(30,-30,0)、B(30,-30,30)位移時程。動力載荷F(t) 的單位為N，表達式為：

(18)

圖5 Lamb表面源算例模型及三參量邊界示意圖

圖6、7分別為A點和B點的位移時程曲線。從圖6、7可以看出，三參量粘彈性邊界的位移時程和兩參數邊界的位移時程基本重合，且均和遠置邊界的近似解吻合得較好，證明了本文使用的三參量粘彈性人工邊界的可行性。

圖6 不同邊界下A點的位移時程曲線

圖7 不同邊界下B點的位移時程曲線

圖8 不同n值下A點的Y向位移時程曲線

1.3.2P波斜入射的波動問題

模型在幾何和邊界上與Lamb表面源問題相同。介質的剪切模量為5.292 GPa，泊松比為0.25，密度為2 700 kg/m3。模擬P波從底部的入射點C(-60,60,-60) 斜入射，入射波與Y軸夾角為60°，入射波與反射波平面和X軸的夾角為45°。入射波在C點引起的位移時程為：

(19)

圖9為觀測點O(0,0,0) 位移時程的解析解和使用粘彈性人工邊界的數值解。如圖9(a)和圖9(b)所示，觀測點O在X、Y方向位移的數值解均和解析解吻合較好，進一步說明了本文使用的三參數粘彈性人工邊界的有效性。

圖9 觀測點O的位移時程

2 埋地輸氣管道有限元模型建立

2.1 管線鋼和土壤的材料參數

管道的材料為X80管線鋼[10]，其彈性模量E為210 GPa，泊松比ν為0.3，密度ρ為7 900 kg/m3。土體使用Drucker-Prager模型，其材料參數見表1[18]。

表1 土體的材料參數

2.2 埋地管道的有限元模型

管道的有限元模型取長度15 m，外徑121.9 cm，壁厚為1.8 cm，埋深(管道中心與地面的距離)為3 m；土體取6 m×6 m×15 m的塊體(Y向為豎直方向；Z向為管道軸向，長度15 m)，其截面的幾何模型如圖10所示。管道的有限元模型如圖11所示，管道和土體均使用Solid186單元，管道與土體的接觸通過定義接觸對實現。模型是在很長的管線中截取的一截，所以在管道軸向的2個邊界面上使用對稱邊界；在模型的底面和2個側面施加三參量粘彈性邊界；模型上表面為自由面。

圖10 埋地管道的幾何模型及監測點示意圖

圖11 埋地管道的有限元模型及邊界

3 地震載荷作用下管道的動態響應

3.1 地震波在模型邊界輸入

埋地管道位于四川江油，根據抗震設計規范，該地區的抗震設防烈度為7°，設計基本地震加速度值為0.15g，設計地震分組為`第二組，建筑場地為Ⅱ類。本文選擇罕遇地震來進行抗震設計，由表2可知，輸入地震加速度的最大值為310 cm/s2。本文使用EI Centro波作為輸入地震波，時長為30 s，時間間隔0.02 s，其加速度時程曲線如圖12所示。地震波的輸入通過式(17)將地震波轉化為等效節點載荷實現。

表2 時程分析時輸入地震加速度的最大值 (cm/s2)

圖12 調整最大加速度后的EI Centro波

3.2 地震波入射角對管道動態響應的影響

由于震源具有很強的不確定性，地震波傳到埋地管道的輸入角也隨之不同。將地震波沿水平和垂直方向分成2個分量，不改變地震波加速度大小，通過調整2個分量的比值，可以得到不同的地震波輸入角度。數值模擬中，輸入角度分別為0°、15°、30°、45°、60°、75°和90°。

選取模型中部截面處管道內壁的底部節點(1號點，節點編號31313)、頂部節點(2號點，節點編號27684)、右端節點(3號點，節點編號27655)和左端節點(4號點，節點編號29513)作為監測點，并分析4個節點在不同入射角下的位移和應力響應。4個節點的具體位置如圖10所示。

3.2.1不同入射角的位移響應

圖13依次給出了入射角為0°、45°、90°時，4個節點的合位移時程曲線。受地震波的作用，4個節點合位移時程曲線的變化趨勢相同，在幅值上有一定的差異。圖13(a)，入射角為0°時，管道底部節點的位移幅值最大；底部的位移幅值最小；圖13(b)，入射角增大到45°時，底部和右端節點、頂部和左端節點的位移幅值分別保持一致，且底部節點和右端節點的位移較大；圖13(c)，入射角為90°時，各節點的位移時程曲線幾乎重合，管道的位移響應呈現出整體一致性。

圖13 不同入射角下各監測點的合位移響應曲線

圖14給出了不同入射角下，管道底部節點的合位移時響應曲線。可以看出，隨著入射角的增大，位移響應也越大；8.62 s，位移達到峰值時，不同入射角的位移差異也達到最大；0°入射角的最大位移為0.61 m，而90°入射角的最大位移達0.131 m。

綜上,管道受地震作用時，管道底部區域的位移響應最大，是管道抗震設計中應著重關注的位置；地震波的入射角越大，管道的位移響應越大，越容易在地震中受損。

圖14 不同入射角下管道底部節點的位移時程曲線

3.2.2不同入射角的應力響應

本文選取Von Mises等效應力分析管道的動態應力響應。選取管道中部截面內壁上的所有節點，并獲取了這些節點在不同入射角地震波作用過程中的最大等效應力。圖15給出了不同輸入角對應的最大等效應力環管道內壁的分布曲線，同時給出了環管道內壁的初始應力分布作為參照。圖15中，極坐標的0°對應圖10中的管道右側的27655號節點；90°對應管道頂部的27684節點；180°對應管道左端的29513號節點；270°對應管道底部的31313號節點。如圖15所示，初始應力呈橢圓形分布，底部和頂部初始應力較大；入射角從0°～45°時，最大等效應力的峰值不斷增大，峰值的位置朝逆時針方向發生偏移；入射角從 45°～90°，最大等效應力的峰值逐漸減小，峰值出現的位置繼續向逆時針方向偏移。

圖15 不同入射角的最大等效應力環向分布圖

為了說明最大等效應力的值和位置隨入射角的變化關系，圖16給出了整個管道在地震過程中的最大等效應力和最大等效應力增量(最大等效應力減去初始時刻的應力)與入射角關系曲線；表3是0°到90°入射角所對應的最大等效應力和最大等效應力增量的位置。

表3 不同入射角管道內最大等效應力和等效應力增量的位置 (°)

如圖16，最大等效應力和應力增量隨地震波入射角的變化趨勢相同，均先增大后減小，最后趨于平穩；最大等效應力和應力增量均在45°時達到最大，值分別為364 MPa和54.4 MPa。如表4，同一入射角的最大等效應力和最大等效應力增量所處的位置基本一致，為同一節點或相鄰節點(管道內壁共有80個節點，相鄰節點位置相差4.5°)；管道內最大等效應力的位置從入射角0°到90°依次為243°、247.5°、247.5°、247.5°、247.5°、265.5°、270°。

圖16 不同入射角下，管道內的最大等效應力及增量曲線

綜上,隨著入射角從0°增大到90°，地震過程中最危險的位置(等效應力最大)由243°(管道底部偏左27°)逐漸轉移到270°(管道正底部)，危險位置的等效應力先增大后減小；地震波0°入射的最大等效應力為332 MPa，對管道強度的威脅最小；地震波45°入射對管道的強度威脅最大，此時的最大等效應力值為364 MPa，位置位于管道內壁247.5°處。

3.3 埋深對管道動態響應的影響

在管道壁厚均為18 mm，地震波輸入角均為90°的情況下，采用管道埋深為3、3.5、4、4.5、5 m研究埋深對管道地震動力響應的影響。由上一節知，地震波90°入射時，管道底部31313號節點的應力和位移響應最大，所以本節著重分析管道底部。

3.3.1不埋深管道的動態位移響應

埋深不同時，31313號節點的Y向位移時程如圖17所示。由圖17知，埋深越大，管道底部的位移響應越小。

圖17 不同埋深管道31313號節點的Y向位移時程曲線

圖18顯示了31313號節點的最大Y向位移與管道埋深的關系。埋深為3 m時，31313號節點的最大位移為0.131 m；埋深增加到5 m時，最大位移減小至0.103 8 m，減小了20.8%。最大位移在相鄰埋深的減小量和減小百分比見表4。埋深從3 m到3.5 m時，管道的位移減小量最大，值為9.02 mm，位移減小的百分比也最大，值為6.88%；隨著埋深繼續增加，位移減小量和減小的百分比均變小；4.5～5 m時，位移減小量僅有4.55 mm，減小位移的效果為3～3.5 m的一半。

圖18 不同埋深的最大Y向位移

表4 相鄰埋深31313號節點的最大位移減小量和減小百分比

由圖17、18和表4可得出，增大管道的埋深可以減小管道的位移響應，增強管道的抗震能力，但位移減小幅度變弱。

3.3.2不同埋深管道的動態應力響應

31313號節點在不同埋深下的等效應力時程曲線如圖19所示。埋深越大，管道的動態應力響應越大，最大等效應力也最大，越容易在地震中被破壞。

為了更詳細地描述埋深對管道應力響應的影響，圖19給出了管道的最大等效應力與埋深的關系，并在表5中給出了相鄰埋深的應力增量和應力增量百分比。由圖20知，隨著管道埋深的增大，最大等效應力呈線性增長，埋深從3 m增至 5 m時，最大等效應力增加了35.5 MPa，平均增長率為2.57%。由表5知，3.5～4 m時，最大應力的增量有最大值10.62 MPa，增加的百分比為3.01%；4.5～5 m時，增大應力的增量有最小值7.79 MPa，增量百分比為2.06%；隨著埋深的增加，管道最大等效應力也平穩增加。

圖19 不同埋深管道31313號節點的等效應力時程曲線

表5 相鄰埋深的最大等效應力量和增量百分比

圖20 不同埋深的最大等效應力曲線

由于管道的應力響應和位移響應隨埋深的變化趨勢相反，所以在進行埋地管道的抗震設計時，需要綜合考慮最大等效應力和最大位移與埋深的關系。在管道強度足夠的前提下，適當增大埋深有利于提升埋地管道的抗震性能。

4 結論

1) 本文建立的三參量粘彈性邊界與傳統兩參量粘彈性邊界具有同等的精度，在描述邊界層的粘彈性時優于傳統的粘彈性邊界。

2) 對于管道內壁的底部、頂部、左右端4個特殊位置，底部的位移響應最大且地震波入射角越大，位移響應最大；對于整個管道，地震波45°入射時，管道的應力響應最大，對管道的強度威脅最大，危險點位于管道內壁247.5°處。

3) 地震波90°入射時，管道埋深越大，管道的位移響應越小，應力響應越大；在管道強度足夠的前提下，適當增大埋深有利于提升埋地管道的抗震性能。

4) 對結構-地基相互作用系統(如埋地輸氣管道)進行動力分析時，邊界是影響計算精度的重要因素。本文提出的三參量粘彈性邊界在具有高精度的同時又能完全描述邊界層的粘彈性，提供的埋地管道地震分析方法和數據對埋地輸氣管道的抗震設計、施工與安全運營具有參考價值。