埋地管線數值模擬中土體計算區域的選取問題研究

劉鵬，陳洪富，孫柏濤.2024.埋地管線數值模擬中土體計算區域的選取問題研究［J］.地震研究，47（1）：018-026，doi：10.20015/j.cnki.ISSN1000-0666.2024.0006.

Liu P，Chen H F，Sun B T.2024.Study on the selection of soil calculation area in numerical simulation of the buried pipeline［J］.Journal of Seismological Research，47（1）：018-026，doi：10.20015/j.cnki.ISSN1000-0666.2024.0006.

摘要：詳細論述了動力無限元基本理論，提出了一種在非一致激勵時無限元邊界面上地震動的輸入方法。在此基礎上，考慮管線-土體之間的黏結滑移效應，采用動力無限元邊界，建立了埋地管線的三維有限元-無限元耦合的精細化模型，對非一致地震波激勵下埋地管線的動力響應進行了數值分析，對管線周圍的土體計算范圍的取值進行了系統研究，給出了管線周圍土體的有效計算寬度W、有效計算長度L和有效計算深度H的取值范圍。結果表明：W的取值與管徑D緊密相關，H可取30～60 m，L應至少包含一個地震動的卓越周期段的波長，此分析結果可滿足大多數工程需求。

關鍵詞：埋地管線；無限元；邊界條件；有限元分析；管線-土體相互作用

中圖分類號：TU311.3"""文獻標識碼：A"""文章編號：1000-0666（2024）01-0018-09

doi：10.20015/j.cnki.ISSN1000-0666.2024.0006

0"引言

埋地管線縱橫交錯、長距離延伸、互相關聯，構成管線與管網。埋地管線在歷次大地震中都有不同程度的震害，甚至會引發次生災害，給國計民生帶來重大損失。由于其隱蔽性和復雜性，管網破壞后在短期內修復困難，嚴重影響震后基礎設施功能恢復的進程。由于埋地管線和管網的抗震試驗受條件限制，其研究范圍和成果有限，采用數值模擬的方法進行分析既能提升效率，又能節約成本，對整套地下管網的地震反應性狀的研究或抗震設計都大有助益。

埋地管線在地震波作用下的反應分析最早由Newmark和Hall（1975）提出，Kennedy等（1977）在Newmark管道模型假定的基礎上進行改進，提出了二維彈性地基梁模型。這2種常用的方法都是將管道周圍的土體對管道的作用簡化為一系列土彈簧的作用，計算精度難以滿足要求。鑒于此，筆者將埋地管線及周圍土體從半無限空間中取出作為一個整體進行研究，但面臨3個需要重點考慮的問題：①地震動如何輸入？②土體的截斷面如何設置人工邊界條件？③管線周圍的土體取值范圍如何確定？

但是這種處理方法計算效率極低，有時甚至低到無法計算的程度。考慮管線-土體相互作用的埋地管線地震反應數值模擬通常只能人為截取足夠大的有限域來模擬實際的半無限空間，此時截斷面必須設置合理的人工邊界，再配合波動輸入方法來消除地震波在截斷面上的反射，諸多學者曾對此開展了深入研究。考慮其計算精度、計算效率以及在通用有限元軟件中的兼容性，應用較為廣泛的人工邊界有黏性邊界、黏彈性邊界、無限元邊界和遠置自由邊界等。其中無限元最早由Ungless（1973）提出，經過Bettess等（1977，1984），Zhang和Zhao（1987）的改進和發展，無限元已經從一維發展到三維，從單向映射擴展到了多向映射。無限元的主要特點為：①局部坐標系中的有限域到整體坐標系中的無限域的映射。即局部坐標趨近于1時，相應整體坐標趨向無窮大，從而實現計算范圍趨于無窮遠點；②無限域上位移衰減過程的描述。即局部坐標趨近于1時，位移趨近零，從而實現無限遠處位移為零的邊界條件，無限元具有較強的魯棒性，其與有限元方法的協調性較好，計算效率和精度都較高。本文利用ABAQUS軟件提供的豐富的無限單元庫（石亦平，周玉蓉，2006；莊茁，2009），在模型截斷面的側邊界和底邊界采用三維動力無限元邊界模擬半無限空間，在非一致地震波激勵下對埋地管線的動力響應進行了數值分析，確定了管線周圍土體的取值范圍。

1"非一致激勵下動力無限元邊界上的地震動輸入方法

1.1"ABAQUS動力無限元基本理論

ABAQUS動力無限元邊界以Lysmer和Kuhlemeyer（1969）的動力響應理論為基礎，由平面體波垂直入射邊界的運動方程進行推導，包含以下基本假設：①有限域邊界區動力反應很小，邊界區介質在線彈性范圍內響應；②無限單元區域介質滿足理想線彈性體假定。則質點運動微分方程可表示為：

根據理想線彈性體假定，由廣義胡克定律可知，質點的應力應變關系可用張量表示為：

σ_ij=λε_kkδ_ij+2Gε_ij" " " " " （2）

現通過一維平面縱波u_x=f₁（x-c_Pt）， u_y=u_z=0垂直入射自由邊界的形式求解分布阻尼系數d_P，如果其在邊界面上完全反射，則應滿足應力為零的自由邊界條件，此時反射波可表示為：

u_x=f₂（x-c_Pt）， u_y=u_z=0" " " " " （9）

根據本節設定的2個基本假設，邊界面上的力和位移均滿足疊加原理，邊界面上的位移、應力和質點振動速度可表示為：

同理可得S波入射時無限元邊界的阻尼系數為：

d_S=ρc_S（13）

1.2"地震動輸入的波動法

從震源出發的地震波，經過在地殼中復雜介質的多次折射、反射，很難確定其入射方向。在整個傳播過程中，地震波的幅值、頻譜組成及其傳播和振動的方向都在不斷改變，且其各個頻率分量的傳播速度也不相同。地震波中所包含的體波、面波也隨到達的時間不同而不斷改變。目前尚難區分在地表記錄到的地震動加速度波形中的各個波型，但場址地表某點的地震動加速度總是可以分解成相互正交的3個分量，通常取為沿地表2個水平方向和垂直于地表的豎直方向。總體上，由于地殼介質的密度由地表向下隨地層深度變深而逐漸增大，根據波在不同介質中傳播的折射和反射定律，由地殼深部向地表傳播的地震波，其入射方向將逐漸變為接近垂直水平地表的豎向方向。對于研究對象為地表以下十幾米到幾十米的地下結構，可以認為地震波在其底面垂直入射。因此，本文地震動的輸入方法均假定地震波從有限域的底邊界垂直入射（胡聿賢，2006；趙武勝等，2010）。

式中：等效節點力的下標為分量方向，上標為節點所在人工邊界的外法線方向，與坐標軸相同為正，反之為負；Δt₁和Δt₂分別為入射P波和地表反射P波的時間延遲；Δt₃和Δt₄為入射S波和地表反射S波的時間延遲。如圖1所示，a為節點到底邊界的距離，L為底邊界到地表的距離。動力計算時可將給出的地震動通過式（15）（16）轉化為節點力，像面荷載一樣以節點力的形式進行施加。

1.3"基于波動法的非一致激勵問題

現階段對空間地震動傳播規律的研究以行波效應為主，基巖面上不同位置的地震激勵是不同步的（圖2）。假設行波傳播方向與管道鋪設方向相同，基底各點受到的地震激勵存在相位差Δt，即地震波從基底一點傳到相鄰的下一點所用的時間為：

式中：Δt₀為最小時間間隔，取值與輸入地震動記錄的時間步長有關。實際輸入的地震動是一系列離散的點，其采樣率一般為Δt₀=0.005或0.01 s，為了保證各點輸入地震動的時間標準與原記錄保持一致，以ΔL_i=cΔt₀為單位對模型在沿著地震動傳播方向上進行等間隔分段，各段ΔL_i內采用一致輸入，相鄰段ΔL_i和ΔL_i+1在時間上相差Δt₀。

2"模型設計

2.1"基本參數

本文研究長地下埋設管線，管線由預制鑄鐵管構成，管材選用X42鑄鐵管，軸向設計拉伸曲線采用Ramberg-Osgood方程（Ramberg，Osgood，1943）擬合，具體參數見表1。三維土體的等效有限域計算尺寸為有效寬度W，土層厚度H，管線分布長度L，土體的本構模型采用線性Drucker-Prager（D-P）模型?？紤]管線-土體接觸區域的應力集中現象，對管線附近的土體離散網格局部加密，三維土體自由場計算模型如圖3所示。

2.2"管線-土體相互作用效應

在管線-土體的相互作用分析中，由于管線和土體剛度存在較大差異，在地震波作用下，管線和土體的交界面上可能發生明顯的相對剪切滑移，為了考慮這種復雜的相互作用效應，在管線-土體交界面上通過設置接觸單元進行模擬。該單元在交界面的法向方向采用硬接觸，即只考慮土體對管道的壓應力，忽略土體對管道的拉應力；在交界面的切向方向采用庫倫摩擦準則計算管線-土體之間的摩擦力為：

τ=μ×P" " " " " （19）

P=WB_bD" " " " " （20）

式中：P為兩接觸面上的法向壓應力；W為土的比重；D為管徑；B_b為管溝頂部寬度；μ為管線-土體之間的摩擦系數，根據《工程地質手冊（第四版）》（常士驃，張蘇民，2007），μ取0.5。當接觸面上的切向剪切應力小于摩擦力時，接觸面沿切向不發生相對滑移；反之，接觸面沿切向發生相對滑移。

2.3"邊界條件

在模型截斷面設置動力無限元邊界，模型的側邊界和底邊界采用三維動力無限元邊界模擬半無限空間，采用CIN3D8單元離散；由于整個模型屬細長結構，根據圣維南原理的局部效應，端部截斷面對管線中間段的受力影響較小，故兩端邊界采用遠置滑動自由邊界，即僅約束界面的y、z方向自由度，不約束x方向自由度。具體設置見表2。

2.4"地震動輸入原則

《油氣輸送管道線路工程抗震技術規范》（GB/T 50470—2017）指出，根據大量震害統計資料，一般場地的地下直埋管道地震動峰值加速度（PGA）大于或等于0.3 g時才開始破壞。為了安全起見，PGA≥0.2 g時，地下直埋管道應進行抗拉伸和抗壓縮驗算。已有研究表明，PGA隨地面下的深度的增加逐漸衰減，但具體的衰減規律尚未達成共識。《蘇聯地震區建筑設計規范》（СНИПⅡ-A 12-69）中規定，地面下100 m深處設計地震加速度可取為地面的50%；《印度標準：結構抗震設計規范》（IS：1893—1984）規定，地面下30 m深處設計地震加速度可減少50%；岡本舜三（1978）建議在地下幾十米深處的設計地震加速度可取為地面的1/2～1/3；《油氣輸送管道線路工程抗震技術規范》（GB/T 50470—2017）建議地下50 m深處的設計地震加速度值取為地面的50%，以上按內插法取用。本文地震動以自由表面的PGA（0.2、0.3、0.4、0.5 g）為基準，y方向輸入時考慮地震波的行波效應，x、z方向一致輸入。加速度最大值按照1（y方向）：0.85（x方向）：0.65（z方向）對其調幅，再按照《油氣輸送管道線路工程抗震技術規范》（GB/T 50470—2017）的建議在基底輸入面折減。

3"埋地管線周圍土體計算范圍取值

管線周圍土體截取足夠大的方盒形有限域、同時在邊界上施加相應的近似約束邊界條件模擬無限域時，對“足夠大”的界定令人困擾：區域較小對數值計算規模的控制很有利，但理論上會帶來較大誤差；區域較大能減小理論誤差，但數值計算規模將以級數倍增加。目前對各類人工邊界計算區域合理取值的討論較少，無限元邊界的位置問題的研究還是空白，以下對埋地管線的有限土體計算區域取值問題進行討論。

3.1"有效計算深度H

場地土層對地震動有放大效應。已有研究表明，埋地管線的破壞率隨著場地覆蓋層厚度的增加而升高，但覆蓋層厚度超出一定范圍后破壞率變化不大。在土層的地震反應分析中，如果以實際基巖面作為地震動輸入面，覆蓋層厚度會變得很大，這在一般工程計算中難以實現。因此本文不在實際基巖面輸入地震波，而是采用在土層的假想基巖面處輸入地震波的方式計算土層的動力反應。工程中常按以下原則確定假想基巖面：①假想基巖面以下各層巖土的剪切波速均不小于500 m/s；②當地面5 m以下存在剪切波速大于其上部各土層剪切波速2.5倍的土層，且該層及下臥各層巖土的剪切波速均不小于400 m/s時，可取該土層頂面為假想基巖面。表3給出了假想基巖面按土層剪切波速和場地覆蓋層厚度劃分的類別。埋地管線所處地質環境多以Ⅱ類場地為主，本文的有效計算深度H取40 m，涵蓋大多數管線所處的地質環境。

3.2"有效計算長度L

L的取值范圍與實際的場址條件有關，其大小至少應反映地震動的主要特性，綜合考慮采用以下方法確定L的最小限值：①求解輸入地震動的Fourier振幅譜的峰值振幅對應的周期T；②求解假想基巖面處土層的剪切波速V_S；③計算該地震動在場地土層中的傳播波長：λ=TV_S；④L應至少包含一個地震動卓越周期段的波長，據此取L≥λ。本文選取EI-Centro地震波，峰值振幅對應的周期為0.42 s，V_S范圍為219.20～489.80 m/s，則L≥205.714 m。

3.3"有效計算寬度W

有效計算寬度W的取值對地震波的反射和散射效應非常敏感，取值過小可能導致巨大誤差，取值過大計算量成倍增加、計算效率較低。為了確定W的合理范圍，本文建立三維精細化模型。為了得到具有規律性的一般結論，對土體模型作如下簡化：土體采用線彈性本構模型；在豎向范圍內不考慮土層的差異性。

實際上，土體進入塑性變形之后對管線的作用力減小，在該假設下由于未考慮塑性變形，計算結果是偏于安全的。本文以埋深1.5 m、管壁厚t=8 mm的埋地管道為例，其土體范圍豎向深度H=40 m、縱向長度L=250 m，分析其有效計算寬度W的取值。選取不同的管徑D，截取不同的土層寬度W，采用動力無限元邊界在模型底邊界輸入4條不同的地震動，將其PGA調幅為0.2、0.3、0.4和0.5 g，比較分析模型中部截面管道頂部的最大軸向應變ε的變化規律。

從不同管徑D、寬度W的管線峰值應變ε變化圖（圖4a）對比發現，當PGA一定時，管道截面的ε隨著D的增大呈明顯的減小趨勢；當D一定時，管道截面的ε隨著PGA的增大呈增大趨勢，并且增幅越來越大，但其曲線的整體走勢基本保持一致；管道截面的ε隨著W的增大逐漸降低直至趨于穩定。當W取值較小時，邊界效應對核心計算區的影響較大，隨著W的逐漸增大，其影響越來越小。這是由于當W取值較小時，無限元邊界不能完全吸收反射波，邊界反射波對計算區域有較大的影響，不可忽略；當W逐漸增大時，反射波對計算區域的影響逐漸減弱。

為了進一步分析W的取值對核心計算區域的影響，本文采用W=40 m（D=0.5，0.7 m）、W=60 m（D=1.0 m）時的管線中部截面某一節點處的軸向峰值應變ε0作為基準值，討論不同PGA取值下管線中部截面的軸向峰值應變ε的相對誤差隨有效計算寬度W的變化情況。其相對誤差可表示為：

式中：εi（i=10～60）為不同有效計算寬度W的管線中部截面同一節點處的軸向峰值應變。

分析圖4b可知：相對誤差隨著W的增大而減小，但不同PGA對計算結果的收斂速度影響不大，最終都能收斂到精確解。當D=0.5 m、W≥20 m時，或者當D=0.7 m、W＞25 m時，相對誤差均可以控制在5%以內；當D=1.0 m、W≥40 m時，相對誤差可以控制在5%以內；當W取值足夠大時，相對誤差就會足夠小。圖4的分析結果表明：W與D有較大關聯，不同的D對W的取值較敏感。當W=40D時，管體的ε明顯趨于穩定，相對誤差均在5%以內，此時可以認為W=40D時可消除邊界反射波對計算區域的影響。

4"結論

本文闡述了無限元基本理論，建立了管線-土體相互作用的三維動力有限元-無限元耦合模型，提出了一種考慮非一致激勵時無限元邊界面上地震動的輸入方法，分析了在地震波作用下埋地管線的動力響應，研究了模型中土體有效計算區域的取值問題，通過多種工況的對比分析，得到了如下結論：

（1）按照土層剪切波速和覆蓋層厚度劃分的假想基巖面作為地震動輸入界面，此時土層的有效深度取值為H=40 m，可涵蓋管線所處的大多數地質環境，具體取值根據管線所處場地的類別決定，一般情況下H=30～60 m。

（2）有效計算長度L的取值與考慮的實際地震波有很大關系，L應至少包含一個地震動的卓越周期段的波長，才可基本反應地震動的時空差異性。

（3）有效計算寬度W與管線的直徑D有很大關聯，W的取值隨著管徑的增大而增大，通過多工況的對比分析，W=40D時可消除邊界反射波對計算區域的影響。

限于文章篇幅，本文以管線周圍土體勻質作為基本假設，未考慮管溝形狀和回填砂土的影響。未來，將考慮管線周圍土體的差異性和管溝形狀的不同對管線地震響應進開展一步研究。

參考文獻：

常士驃，張蘇民.2007.工程地質手冊（第四版）［M］.北京：中國建筑工業出版社.Chang S B，Zhang S M.2007.Engineering geology handbook（4th Edition）［M］.Beijing：China Architecture amp; Building Press.（in Chinese）

杜修力，趙密，王進廷.2006.近場波動模擬的人工應力邊界條件［J］.力學學報，38（1）：49-56.Du X L，Zhao M，Wang J T.2006.A stress artificial boundary in FEA for near-field wave problem［J］.Chinese Journal of Theoretical and Applied Mechanics，38（1）：49-56.（in Chinese）

杜修力.2009.工程波動理論與方法［M］.北京：科學出版社.Du X L.2009.Theories and methods of wave motion for engineering［M］.Beijing：Science Press.（in Chinese）

岡本舜三.1978.抗震工程學［M］.北京：中國建筑工業出版社.Okamoto S.1978.Aseismic engineering［M］.Beijing：China Architecture amp; Building Press.（in Chinese）

谷音，劉晶波，杜義欣.2007.三維一致黏彈性人工邊界及黏彈性邊界單元［J］.工程力學，24（12）：31-37.Gu Y，Liu J B，Du Y X.2007.3D Consistent viscous-spring artificial boundary and viscous-spring boundary element.［J］.Engineering Mechanics，24（12）：31-37.（in Chinese）

胡聿賢.2006.地震工程學（第2版）［M］.北京：地震出版社.Hu Y X.2006.Earthquake engineering（2nd Edition）［M］.Beijing：Seismological Press.（in Chinese）

劉晶波，古音，杜義欣.2006.一致黏彈性人工邊界及黏彈性邊界單元［J］.巖土工程學報，28（9）：1070-1075.Liu J B，Gu Y，Du Y X.2006.Consistent viscous-spring artificial boundary and viscous-spring boundary element［J］.Chinese Journal of Geotechnical Engineering，28（9）：1070-1075.（in Chinese）

劉晶波，王振宇，杜修力，等.2005.波動問題中的三維時域黏彈性人工邊界［J］.工程力學，22（6）：46-51.Liu J B，Wang Z Y，Du X L.2005.Three-dimensional visco-elastic artificial boundaries in time domain for wave motion problems［J］.Engineering Mechanics，22（6）：46-51.（in Chinese）

石亦平，周玉蓉.2006.ABAQUS有限元分析實例詳解［M］北京：機械工業出版社.Shi Y P，Zhou Y R.2006.Detailed explanation of ABAQUS finite element analysis example［M］.Beijing：China Machine Press.（in Chinese）

趙建峰，杜修力，韓強，等.2007.外源波動問題數值模擬的一種實現方式［J］.工程力學，24（4）：52-58.Zhao J F，Du X L，Han Q，et al.2007.An approach to numerical simulation for external source wave motion［J］.Engineering Mechanics，24（4）：52-58.（in Chinese）

趙武勝，陳衛忠，黃勝，等.2010.地下結構抗震分析中地震動輸入方法研究［J］.土木工程學報.43（S1）：554-559.Zhao W S，Chen W Z，Huang S，et al.2010.Research on method of seismic motion input in aseismic analysis for underground structure［J］.China Civil Engineering Journal，43（S1）：554-559.（in Chinese）

莊茁.2009.基于ABAQUS的有限元分析和應用［M］北京：清華大學出版社.Zhuang Z.2009.Finite element analysis and application based on ABAQUS［M］.Beijing：Tsinghua University Press.（in Chinese）

CHИПⅡ-A 12-69，蘇聯地震區建筑設計規范［S］.

CHИПⅡ-A 12-69，Building design code for earthquake area in USSR［S］.（in Chinese）

GB/T 50470—2017，油氣輸送管道線路工程抗震技術規范［S］.

GB/T 50470—2017，Seismic technical code for oil and gas transmission pipeline engineering［S］.（in Chinese）

GB 50011—2010，建筑抗震設計規范［S］.

GB 50011—2010，Code for seismic design of buildings［S］.（in Chinese）

IS：1893—1984（2003），印度標準：結構抗震設計規范［S］.

IS：1893—1984（2003），Indian standard：Criteria for earthquake resistant design of structures［S］.（in Chinese）

Bettess P，Eemson C，Chiam T C.1984.A new mapped infinite element for exterior wave problems［M］.New York：John Wiley amp; Sons.

Bettess P，Zienkiewicz.1977.Diffraction and refraction of surface waves using finite and infinite elements［J］.International Journal for Numerical Methods in Engineering，11（6）：1271-1290.

Kennedy R P，Chow A W，William R A.1977.Fault movement effects on buried oil pipeline［J］.Journal of Transportation Engineering，ASCE，103（6）：617-633.

Lysmer J，Kuhlemeyer R L.1969.Finite dynamic model for infinite media［J］.Journal of Engineering Mechanics，ASCE，95（4）：859-877.

Newmark N M，Hall W J.1975.Pipeline design to resist large fault displacement［C］//Earthquake Engineering Research Institute.Proceedings of US National Conference on Earthquake Engineering，416-425.

Ramberg W，Osgood W R.1943.Description of stress-strain curves by three parameters［R］.National Advisory Committee For Aeronautics，Washington D C，Technical Note No.902.

Ungless R F.1973.An infinite finite element［D］.Vancouver：the University of British Columbia.

Zhang C H，Zhao C B.1987.Coupling method of finite and infinite elements for strip foundation wave problems［J］.Earthquake Engineering and Structural Dynamics，15（7）：839-851.

Study on the Selection of Soil Calculation Area in Numerical Simulationof the Buried Pipeline

LIU Peng^1，2，CHEN Hongfu^1，2，SUN Baitao^1，2

（1.Key Laboratory of Earthquake Engineering and Engineering Vibration，Institute of Engineering Mechanics，China Earthquake Administration，Harbin 150080，Heilongjiang，China）

（2.Key Laboratory of Earthquake Disaster Mitigation，Ministry of Emergency Management，Harbin 150080，Heilongjiang，China）

Abstract

In this paper，the basic theory of dynamic infinite element is discussed，and an input method for ground motion at the boundary surface of the infinite element by non-uniform excitation is proposed.Then，on condition that the bond slip effect between pipelines and soil is involved，a refined model of buried pipelines，which couples three-dimensional finite element with infinite element，is established by using dynamic infinite element boundaries.The dynamic response of buried pipelines by the non-uniform seismic wave excitation is numerically analyzed.The values of the calculation range of the soil around the pipeline are systematically studied，and the effective width W，the effective length L，and the effective depth H of the soil around the pipeline are given.The results indicate that the value of W is closely related to the pipe diameter D；H can be taken between 30-60 m，and L should include at least one wavelength of the predominant period of the ground motion.In this case，the results can meet most engineering requirements.

Keywords：buried pipeline；infinite element；boundary condition；finite element analysis；pipe-soil interaction

*收稿日期：2023-03-17.

基金項目：黑龍江省重點研發計劃項目（GA22C001）；國家重點研發計劃課題（2019YFC1509301）.

第一作者簡介：劉"鵬（1989-），博士研究生在讀，主要從事地下管網抗震研究.E-mail：liupeng265@buaa.edu.cn.

?通信作者簡介：陳洪富（1983-），副研究員，博士，主要從事地震災害風險評估和抗震加固研究.E-mail：chen-hongfu@163.com.