ABAQUS技術在地下管道動力分析中的應用

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

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

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

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

1 常用人工邊界

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

1.1 黏性人工邊界

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

(1)

(2)

1.2 黏彈性人工邊界

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

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

在土-地下結構動力分析中推導出人工邊界條件，假定散射波為柱面波，柱面波運動方程為[5-6]：

(3)

(4)

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

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

(5)

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

(6)

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

(7)

Db=ρ·cs

(8)

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

1.3 固定人工邊界

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

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

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

(9)

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

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

2 相互作用

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

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

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

圖3 罰函數的接觸壓力與接觸間隙關系Fig.3 Relationship of contact pressure and gap of penalty function

在ABAQUS軟件中，計算極限剪應力[7]采用公式為：

τcr=μp

(10)

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

(11)

式中：μ為摩擦系數，σy為土壤的 mises 屈服應力。

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

3 實例

3.1 有限元計算模型

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

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

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

圖4 有限元模型網格劃分Fig.4 Grid meshing of finite element model

3.2 計算方法

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

表2 動力反應分析流程Tab.2 Process of dynamic response analysis

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

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

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

3.3 邊界與地震載荷

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

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

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

地震荷載采用人工地震動合成波。為簡化計算，地震波選用沈陽市地鐵四號線文貿路站-文東街站區間的位移時程荷載，彈性波垂直從下方入射，直接作用在模型底部截面上，加載地震波時位移時程選用50年10%超越概率(設防烈度)時程曲線(見圖6)，地震持續時間取為90 s，當地地震設防烈度為7度，設防地震加速度最大值不超過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 計算結果分析

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

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

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

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

表4 最大等效應力分布情況Tab.4 Distribution of maximum equivalent stress

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

表5 四種工況計算結果對比Tab.5 Comparison of the calculation results under four working conditions

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

5 結論

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