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

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

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

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

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

1 三參量粘彈性人工邊界

1.1 現(xiàn)有的粘彈性人工邊界及其缺點

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

(1)

(2)

圖1 管道結(jié)構(gòu)的兩參數(shù)粘彈性邊界地震模型示意圖

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

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

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

1.2.1有限地震模型法向邊界等效參數(shù)的推導(dǎo)

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

(3)

(4)

式中:U為邊界節(jié)點的位移；σ為位移U在人工邊界處產(chǎn)生的應(yīng)力；K1和K2分別為圖2中對應(yīng)位置彈簧的剛度；U1為與邊界節(jié)點相鄰彈簧的位移；U2為遠離邊界節(jié)點彈簧的位移；UC為阻尼器的位移；UM為集中質(zhì)量點的位移；C為阻尼器的阻尼系數(shù)；M為集中質(zhì)量。

圖2 邊界法向的三參數(shù)粘彈性物理系統(tǒng)示意圖

由式(3)和(4)可得邊界節(jié)點應(yīng)力和位移的關(guān)系：

(5)

(6)

球坐標(biāo)系中，P波波陣面上法向應(yīng)力與位移的關(guān)系[16]為：

(7)

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

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

(8)

1.2.2有限地震模型切向邊界等效參數(shù)的推導(dǎo)

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

圖3 邊界切向的三參數(shù)粘彈性物理系統(tǒng)

Kelvin-Voigt模型應(yīng)力與位移的關(guān)系為：

(9)

(10)

球坐標(biāo)系中，S波波陣面上切向應(yīng)力與位移的關(guān)系[16]為：

(11)

式中：cS為介質(zhì)的S波波速。

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

(12)

1.2.3數(shù)值模擬技術(shù)中的參數(shù)設(shè)置

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

(13)

(14)

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

數(shù)值模擬時，將地震波加速度轉(zhuǎn)換為人工粘彈性邊界上的等效節(jié)點載荷來實現(xiàn)地震波的輸入。在邊界上應(yīng)力和位移滿足圖4所示的關(guān)系：

(15)

(16)

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

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

(17)

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

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

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

1.3.1三維Lamb表面源問題

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

(18)

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

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

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

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

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

1.3.2P波斜入射的波動問題

模型在幾何和邊界上與Lamb表面源問題相同。介質(zhì)的剪切模量為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) 位移時程的解析解和使用粘彈性人工邊界的數(shù)值解。如圖9(a)和圖9(b)所示，觀測點O在X、Y方向位移的數(shù)值解均和解析解吻合較好，進一步說明了本文使用的三參數(shù)粘彈性人工邊界的有效性。

圖9 觀測點O的位移時程

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

2.1 管線鋼和土壤的材料參數(shù)

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

表1 土體的材料參數(shù)

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單元，管道與土體的接觸通過定義接觸對實現(xiàn)。模型是在很長的管線中截取的一截，所以在管道軸向的2個邊界面上使用對稱邊界；在模型的底面和2個側(cè)面施加三參量粘彈性邊界；模型上表面為自由面。

圖10 埋地管道的幾何模型及監(jiān)測點示意圖

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

3 地震載荷作用下管道的動態(tài)響應(yīng)

3.1 地震波在模型邊界輸入

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

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

圖12 調(diào)整最大加速度后的EI Centro波

3.2 地震波入射角對管道動態(tài)響應(yīng)的影響

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

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

3.2.1不同入射角的位移響應(yīng)

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

圖13 不同入射角下各監(jiān)測點的合位移響應(yīng)曲線

圖14給出了不同入射角下，管道底部節(jié)點的合位移時響應(yīng)曲線?？梢钥闯?，隨著入射角的增大，位移響應(yīng)也越大；8.62 s，位移達到峰值時，不同入射角的位移差異也達到最大；0°入射角的最大位移為0.61 m，而90°入射角的最大位移達0.131 m。

綜上,管道受地震作用時，管道底部區(qū)域的位移響應(yīng)最大，是管道抗震設(shè)計中應(yīng)著重關(guān)注的位置；地震波的入射角越大，管道的位移響應(yīng)越大，越容易在地震中受損。

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

3.2.2不同入射角的應(yīng)力響應(yīng)

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

圖15 不同入射角的最大等效應(yīng)力環(huán)向分布圖

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

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

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

圖16 不同入射角下，管道內(nèi)的最大等效應(yīng)力及增量曲線

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

3.3 埋深對管道動態(tài)響應(yīng)的影響

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

3.3.1不埋深管道的動態(tài)位移響應(yīng)

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

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

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

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

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

由圖17、18和表4可得出，增大管道的埋深可以減小管道的位移響應(yīng)，增強管道的抗震能力，但位移減小幅度變?nèi)酢?/p>

3.3.2不同埋深管道的動態(tài)應(yīng)力響應(yīng)

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

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

圖19 不同埋深管道31313號節(jié)點的等效應(yīng)力時程曲線

表5 相鄰埋深的最大等效應(yīng)力量和增量百分比

圖20 不同埋深的最大等效應(yīng)力曲線

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

4 結(jié)論

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

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

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

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