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

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

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

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

1 常用人工邊界

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

1.1 黏性人工邊界

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

(1)

(2)

1.2 黏彈性人工邊界

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

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

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

(3)

(4)

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

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

(5)

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

(6)

方程(6)的物理意義是：如果在半徑rb處截斷介質(zhì)設(shè)定人工邊界，則需要在截斷的人工邊界上施加一定剛性系數(shù)的線性彈簧和一定阻尼系數(shù)的黏性阻尼器兩類物理元件，將方程(6)中的系數(shù)用彈簧剛度系數(shù)Kb和阻尼器的阻尼系數(shù)Db來代替，這樣就得到了黏彈性邊界的邊界條件，即為：

(7)

Db=ρ·cs

(8)

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

1.3 固定人工邊界

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

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

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

(9)

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

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

2 相互作用

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

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

τcr=μp

(10)

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

(11)

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

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

3 實例

3.1 有限元計算模型

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

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

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

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

3.2 計算方法

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

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

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

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

管土相互作用采用接觸面法，接觸類型為surface-to-surface contact(standard)，滑動方式采用有限滑移。接觸屬性中法向接觸采用罰函數(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個方向的線位移和3個方向的角位移；導(dǎo)向架-限制管道2個方向的線位移；支托架-限制管道1個方向的線位移。

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

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

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

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

4 計算結(jié)果分析

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

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

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

通過地應(yīng)力平衡分析獲得動荷載加載前的土體應(yīng)力云圖(見圖8)。土體在靜力分析時的Mises應(yīng)力主要集中在下部，應(yīng)力云圖一層層由大到小從下方向上分布，其值在0.05 MPa左右。同時，從圖7上可以看出，管道在本次人工地震作用下出現(xiàn)最大應(yīng)力的時間在step time在4～10 s。最大等效應(yīng)力的單元、節(jié)點分布情況見表4，四種工況計算結(jié)果對比情況見表5，其對應(yīng)時間的管道等效應(yīng)力云圖見圖8。管道應(yīng)力云圖顯示在靠近固定端部的應(yīng)力較大，遠離固定端部的較小，最大值在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 四種工況計算結(jié)果對比Tab.5 Comparison of the calculation results under four working conditions

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

5 結(jié)論

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