國 艷，張文靜，肖 遙，李林侗

(遼寧省地震局，沈陽 110034)

地下工程的快速發(fā)展給人們帶來了便利，但其抗震安全性也越發(fā)重要。地下結(jié)構(gòu)建設(shè)花費(fèi)高，修復(fù)難，一旦破壞，將會(huì)產(chǎn)生巨大的經(jīng)濟(jì)損失與社會(huì)影響。相關(guān)研究指出[1]，由于管道自身重量較小，管道周圍土體對(duì)其約束較大，其應(yīng)力主要由周圍土體的相對(duì)位移引起，在分析地下管道時(shí)，對(duì)于橫截面積不大的直管，彎曲應(yīng)力較小，可僅考慮軸向應(yīng)變。在進(jìn)行地下結(jié)構(gòu)抗震性能分析時(shí)，可利用ABAQUS有限元分析軟件。對(duì)計(jì)算結(jié)果的準(zhǔn)確性和收斂性起關(guān)鍵作用的因素包括荷載與邊界條件、相互作用、分析方法等。

為簡化分析模型，假設(shè)土體為各項(xiàng)同性，靜力分析時(shí)僅考慮自重和固定人工邊界，動(dòng)力分析時(shí)改變荷載和邊界，采用黏彈性邊界和地震荷載，在土體與結(jié)構(gòu)間相互作用時(shí)采用接觸面法來進(jìn)行模擬。

1 常用人工邊界

應(yīng)用數(shù)值方法研究地震波在介質(zhì)中的傳播，考慮土-地下結(jié)構(gòu)動(dòng)力相互作用(SSI)，許多學(xué)者提出多種人工邊界理論。進(jìn)行數(shù)值模擬時(shí)，人工邊界方法非常多，包括靜力人工邊界、動(dòng)力人工邊界、固定邊界、滾軸邊界、黏性人工邊界、黏彈性人工邊界、一致邊界、旁軸邊界、透射邊界、自由邊界、剛性邊界、吸收邊界、顯式人工邊界、整體人工邊界、局部人工邊界、離散人工邊界、隱式人工邊界等。

1.1 黏性人工邊界

黏性人工邊界理論最早是1969年John Lyrsmer[2]提出來的，主要是研究無限域的動(dòng)力分析數(shù)值方法，適用于瞬態(tài)和穩(wěn)態(tài)問題。將無限域替換為有限域，該區(qū)域受邊界條件的約束，邊界條件可模擬能量吸收邊界。通過有限元法進(jìn)行分析，在人工邊界設(shè)置阻尼器用以吸收能量，解決應(yīng)用于地基振動(dòng)問題。黏性人工邊界條件可以表示為：

(1)

(2)

1.2 黏彈性人工邊界

黏彈性人工邊界理論是1994年Deeks[3]提出來的，采用的推導(dǎo)過程與黏性邊界類似，實(shí)質(zhì)是在人工邊界上布置具有一定阻尼系數(shù)的阻尼器和一定剛性系數(shù)的線性彈簧[4]，示意圖見圖1。

圖1 黏彈性人工邊界示意圖Fig.1 Sketch of viscoelastic artificial boundary

在土-地下結(jié)構(gòu)動(dòng)力分析中推導(dǎo)出人工邊界條件，假定散射波為柱面波，柱面波運(yùn)動(dòng)方程為[5-6]：

(3)

(4)

式中：cs為剪切波速；G為剪切模量；ρ為介質(zhì)密度。

柱面波可以近似的表示為：

(5)

利用方程(3)(5)進(jìn)一步推導(dǎo)出任一半徑rb處的微元面上剪應(yīng)力同該處速度和位移關(guān)系可以表示為如下形式：

(6)

方程(6)的物理意義是：如果在半徑rb處截?cái)嘟橘|(zhì)設(shè)定人工邊界，則需要在截?cái)嗟娜斯み吔缟鲜┘右欢▌傂韵禂?shù)的線性彈簧和一定阻尼系數(shù)的黏性阻尼器兩類物理元件，將方程(6)中的系數(shù)用彈簧剛度系數(shù)Kb和阻尼器的阻尼系數(shù)Db來代替，這樣就得到了黏彈性邊界的邊界條件，即為：

(7)

Db=ρ·cs

(8)

其中：rb是黏彈性人工邊界上節(jié)點(diǎn)至波源的距離；Kb和Db分別是人工邊界上的彈簧剛度系數(shù)和阻尼器的阻尼系數(shù)。如果r無窮大，則黏彈性邊界退化為黏性邊界。

1.3 固定人工邊界

固定邊界條件也稱Dirichlet邊界條件/狄利克雷邊界條件/第一類邊界條件，是假定在人工邊界上位移和轉(zhuǎn)角都為零，那么邊界上能量沒有被吸收全部反射，反射系數(shù)為-1，給出未知函數(shù)在人工邊界上的數(shù)值。人工邊界需要土體一般為地下結(jié)構(gòu)的20～30倍，參與有限元計(jì)算的土體范圍相比無限元邊界和遠(yuǎn)置邊界，有限元離散單元數(shù)少很多，計(jì)算工作量小很多。

對(duì)于偏微分方程，?2φ+φ=0，域Ω?Rn的狄利克雷邊界條件采取形式：

φ(x)=f(x)，?x∈?Ω

(9)

其中f(x)是在邊界?Ω中定義的已知函數(shù)。

Dirichlet邊界條件常用于機(jī)械工程和土木工程(梁理論)中，梁的一端保持在空間中的固定位置。在多物理場仿真ABAQUS軟件中，Symmetry/Antisymmetry/Encastre(對(duì)稱/非對(duì)稱/端部固定邊界條件)、displacement/rotation(位移和旋轉(zhuǎn)邊界條件)取值為零，就可以固定相應(yīng)的位移或轉(zhuǎn)角自由度，達(dá)到邊界固定的目的。

2 相互作用

在進(jìn)行有限元分析過程中，對(duì)結(jié)構(gòu)與土體間相互作用的數(shù)值模擬，通常根據(jù)所掌握的數(shù)據(jù)資料采用土彈簧法、PSI單元法、接觸面法等[7]。在ABAQUS軟件的接觸數(shù)值模擬中采用接觸面法，建立管道和土體的實(shí)體有限元模型，在各個(gè)構(gòu)件上建立接觸表面集(surface集)，在各個(gè)表面集上利用interaction模塊來定義管-土接觸對(duì)、接觸對(duì)屬性，采用單純的主-從接觸算法，最后施加地震時(shí)程載荷，完成數(shù)值分析計(jì)算。

在ABAQUS軟件中，結(jié)構(gòu)與土體接觸面之間的相互作用包括法向行為和切向行為[8]，當(dāng)兩個(gè)構(gòu)件間接觸壓力為零或負(fù)值時(shí)，接觸面分離，無法傳遞法向壓力，所以約束消失，這時(shí)法向行為的接觸稱為“硬”接觸(圖2)，而兩個(gè)接觸面間的相對(duì)運(yùn)動(dòng)(滑動(dòng))和摩擦剪應(yīng)力則稱為相互作用的切向行為。通常有限元數(shù)值模擬真實(shí)的摩擦行為十分復(fù)雜，簡化分析后采用默認(rèn)的罰摩擦公式(圖3)，對(duì)于理想的粘結(jié)-滑動(dòng)摩擦行為，使用拉格朗日摩擦公式，產(chǎn)生相對(duì)滑動(dòng)時(shí)需要考慮靜摩擦系數(shù)和動(dòng)摩擦系數(shù)，這時(shí)可以采用指數(shù)衰減關(guān)系來模擬。

圖2 硬接觸的接觸壓力與接觸間隙關(guān)系Fig.2 Relationship of contact pressure and gap of hard contact

圖3 罰函數(shù)的接觸壓力與接觸間隙關(guān)系Fig.3 Relationship of contact pressure and gap of penalty function

在ABAQUS軟件中，計(jì)算極限剪應(yīng)力[7]采用公式為：

τcr=μp

(10)

計(jì)算時(shí)，需要指定一個(gè)所允許的最大剪應(yīng)力[7]采用以下公式：

(11)

式中：μ為摩擦系數(shù)，σy為土壤的 mises 屈服應(yīng)力。

主從面選取應(yīng)遵循：主控表面應(yīng)為剛度較大的材料組成，主控表面網(wǎng)格劃分比較稀疏?；谝陨纤枳裱瓌t，在計(jì)算分析過程中，剛度較大的管道外表面作為主接觸面，剛度相對(duì)較小的土體接觸管道部分作為從接觸面，形成一對(duì)接觸。

3 實(shí)例

3.1 有限元計(jì)算模型

進(jìn)行有限元分析計(jì)算通常采用有限元軟件ABAQUS，該軟件能夠求解分析各類結(jié)構(gòu)的靜力、動(dòng)力、線性、非線性問題。鋼管采用各項(xiàng)同性的彈塑性模型，采用Mises屈服準(zhǔn)則，通常工程中鋼管管道長度3～12 m。本研究中的管道長度為6 m，管道與土體計(jì)算模型參數(shù)見表1。

表1 管道與土體模型參數(shù)Tab.1 Model parameters of pipe and soil

基于管道模型、土體模型、管土之間相互作用建立有限元模型，進(jìn)行網(wǎng)格劃分。有限元網(wǎng)格劃分時(shí)，與管道接觸的土體進(jìn)行網(wǎng)格加密。為了避免主從面穿透問題，劃分網(wǎng)格時(shí)在管道接觸的土體表面網(wǎng)格數(shù)量與管道表面網(wǎng)格數(shù)量相同。網(wǎng)格控制采用中軸算法、掃掠技術(shù)劃分。為了縮短計(jì)算時(shí)間，有限計(jì)算區(qū)域內(nèi)，采用細(xì)網(wǎng)格劃分的線性縮減積分單元。模型中管道與土體計(jì)算單元的單元族均采用3D stress類型， 六面體八節(jié)點(diǎn)等參元進(jìn)行離散，管道單元類型采用實(shí)體模型C3D8R，土體單元類型采用實(shí)體模型C3D8類型。模型總體單元數(shù)量17 200個(gè)，總體節(jié)點(diǎn)數(shù)量19 803個(gè)。其中，土體模型單元數(shù)量16 880個(gè)、節(jié)點(diǎn)數(shù)量19 131個(gè)，管道模型單元數(shù)量320個(gè)、節(jié)點(diǎn)數(shù)量672個(gè)(如圖4所示)。

圖4 有限元模型網(wǎng)格劃分Fig.4 Grid meshing of finite element model

3.2 計(jì)算方法

在ABAQUS中，通常的動(dòng)力反應(yīng)分析方法有動(dòng)態(tài)顯示分析方法(Dynamic Explicit)和動(dòng)態(tài)隱式分析方法(Dynamic Implicit)。常用的分析流程可歸納為3種(見表2)。

表2 動(dòng)力反應(yīng)分析流程Tab.2 Process of dynamic response analysis

計(jì)算方法采用上述分析流程1提供的方法。對(duì)管道采用靜力方法分析，計(jì)算區(qū)域離散單元的尺寸為0.5，靜力分析時(shí)間步長1 s，獲得管道在自重作用下的應(yīng)力場；采用動(dòng)力方法分析，施加地震位移時(shí)程荷載，管道動(dòng)力分析持續(xù)時(shí)間為50 s，計(jì)算管道動(dòng)力響應(yīng)。

整個(gè)載荷歷程劃分為2個(gè)分析步，在第一個(gè)分析步(Static，General)中僅施加靜態(tài)恒載荷，即自重載荷來計(jì)算靜力分析；在第二個(gè)分析步(Dynamic，Implicit)中施加XYZ三個(gè)方向的地震時(shí)程荷載計(jì)算管道的動(dòng)力響應(yīng)。

管土相互作用采用接觸面法，接觸類型為surface-to-surface contact(standard)，滑動(dòng)方式采用有限滑移。接觸屬性中法向接觸采用罰函數(shù)方法，摩擦系數(shù)為0.3，切向接觸采用硬接觸，阻尼方式在linear(standard only)中取值為damp=0.02，clearance=0.0；damp=0.0，clearance=0.01。

3.3 邊界與地震載荷

土體介質(zhì)彈性模量E=20 MPa，密度ρ=2 000 kg/m3，泊松比v=0.25，那么橫波波速Cs=63.2 m/s，縱波波速Cp=109.5 m/s。通常管道固定方式分為3種：固定架-限制管道3個(gè)方向的線位移和3個(gè)方向的角位移；導(dǎo)向架-限制管道2個(gè)方向的線位移；支托架-限制管道1個(gè)方向的線位移。

土體人工邊界靜力分析時(shí)采用固定邊界，水平方向邊界u1=u3=0，豎直方向邊界u2=0。動(dòng)力分析時(shí)采用黏彈性邊界，在人工邊界網(wǎng)格節(jié)點(diǎn)切線和法線上施加阻尼器和彈簧。本研究選取阻尼器和彈簧類型為連接點(diǎn)與地面方式(見圖5)。

圖5 有限元模型邊界Fig.5 Boundary of finite element model

地震荷載采用人工地震動(dòng)合成波。為簡化計(jì)算，地震波選用沈陽市地鐵四號(hào)線文貿(mào)路站-文東街站區(qū)間的位移時(shí)程荷載，彈性波垂直從下方入射，直接作用在模型底部截面上，加載地震波時(shí)位移時(shí)程選用50年10%超越概率(設(shè)防烈度)時(shí)程曲線(見圖6)，地震持續(xù)時(shí)間取為90 s，當(dāng)?shù)氐卣鹪O(shè)防烈度為7度，設(shè)防地震加速度最大值不超過100 cm/s2。本次計(jì)算時(shí)程荷載加速度最大值X方向52.4 cm/s2、Y方向53.1 cm/s2、Z方向53.6 cm/s2。

圖6 地震荷載時(shí)程曲線Fig.6 Time-history curves of earthquake load

4 計(jì)算結(jié)果分析

為了分析管道不同固定方式的動(dòng)力響應(yīng)，分別對(duì)4種工況(見表3)下管道的受力性能進(jìn)行分析，且由于管道自重相對(duì)較小，管道周圍土體對(duì)其約束較大，其應(yīng)力主要由周圍土體的相對(duì)位移引起。對(duì)于橫截面積不大的直管進(jìn)行有限元?jiǎng)恿Ψ磻?yīng)分析，在計(jì)算的應(yīng)力結(jié)果中找到整個(gè)時(shí)間歷程中應(yīng)力最大的單元，提取該單元的整個(gè)時(shí)間歷程曲線(見圖7)。

表3 工況列表Tab.3 List of working conditions

圖7 四種工況最大等效應(yīng)力曲線圖Fig.7 Maximum equivalent stress curve under four working conditions

通過地應(yīng)力平衡分析獲得動(dòng)荷載加載前的土體應(yīng)力云圖(見圖8)。土體在靜力分析時(shí)的Mises應(yīng)力主要集中在下部，應(yīng)力云圖一層層由大到小從下方向上分布，其值在0.05 MPa左右。同時(shí)，從圖7上可以看出，管道在本次人工地震作用下出現(xiàn)最大應(yīng)力的時(shí)間在step time在4～10 s。最大等效應(yīng)力的單元、節(jié)點(diǎn)分布情況見表4，四種工況計(jì)算結(jié)果對(duì)比情況見表5，其對(duì)應(yīng)時(shí)間的管道等效應(yīng)力云圖見圖8。管道應(yīng)力云圖顯示在靠近固定端部的應(yīng)力較大，遠(yuǎn)離固定端部的較小，最大值在140～420 MPa。

表4 最大等效應(yīng)力分布情況Tab.4 Distribution of maximum equivalent stress

圖8 有限元模型受力分析- Mises應(yīng)力Fig.8 Force analysis of finite element model-Mises stress

表5 四種工況計(jì)算結(jié)果對(duì)比Tab.5 Comparison of the calculation results under four working conditions

當(dāng)采用實(shí)體模型分析管道動(dòng)力響應(yīng)時(shí)，模型計(jì)算得到了相對(duì)穩(wěn)定、易收斂的結(jié)果。計(jì)算結(jié)果對(duì)比顯示，直埋無固定管道模型響應(yīng)最小，最大節(jié)點(diǎn)應(yīng)力僅為屈服應(yīng)力的35.41%，導(dǎo)向架管道模型響應(yīng)最大，最大節(jié)點(diǎn)應(yīng)力為屈服應(yīng)力99%。鑒于地震安全性和工程成本造價(jià)考慮，管道鋪設(shè)最好的方式是直埋無固定管道，其次是支托架管道、固定架管道，最后是導(dǎo)向架管道。實(shí)際工程中，可參考當(dāng)?shù)氐牡刭|(zhì)情況酌情選取支托架、固定架、導(dǎo)向架。通過有限元方法對(duì)埋地管道進(jìn)行應(yīng)力分析結(jié)果對(duì)比，選取適合的管道固定方式，既能減少工程成本造價(jià)，也能保障埋地管道抗震設(shè)防能力，為實(shí)際工程提供相應(yīng)的技術(shù)佐證。

5 結(jié)論

通過運(yùn)用大型通用軟件ABAQUS，對(duì)幾種固定方式管道的有限元數(shù)值模擬結(jié)果進(jìn)行比對(duì)，結(jié)果表明，用ABAQUS計(jì)算管土地下結(jié)構(gòu)的動(dòng)力分析是可行且符合實(shí)際的。通過以上算例可以看出，通過進(jìn)行管道動(dòng)力響應(yīng)分析，可以查找出最有利的管道固定方式，同時(shí)還能找出模型最大應(yīng)力點(diǎn)，也就是管道的薄弱位置，為城市輔助決策的制定提供有利依據(jù)。