二維黏彈性人工邊界單元及地震波輸入在ANSYS中的實現

尹訓強,羅 勇,王桂萱

( 大連大學 土木工程技術研究與開發中心，遼寧 大連 116622 )

對無限域或半無限域進行地震響應分析時，往往采用有限元法。該方法通常采取的處理措施是在無限域或半無限域空間中人為地截取一定范圍的有限域，進而可以將無限區域轉化為有限區域，然后在有限域的基礎上，對其進行離散法處理，進而無限域或半無限域最終被轉化為有限個單元。同時為了達到近似處理的效果，通常在截取的有限域的邊界上人為的施加約束，但是這一處理的缺點是在處理過程中，不僅沒有考慮遠場無限地基的輻射阻尼效應，同時半空間無限地基彈性恢復性能也被忽略。由于截取的有限域的范圍存在不確定性，隨著有限域截取范圍的增加，自由度也隨之增加。基于此離散范圍也更加寬廣，進而網格劃分的單元數及結點數成倍增加，因此，在地震動的動力計算中工作量也進一步被增加，對計算機的要求也更高，進而，在動力學中非線性問題的求解難度進一步增大。目前，針對該關鍵性問題處理的諸多措施中，較好的處理方法是在截取的有限域邊界處設置局部人工邊界。目前在工程中較常用的局部邊界主要有黏性邊界[1]、透射邊界[2]、黏彈性邊界[3]等。黏性邊界的優點是物理概念清晰，在實際動力計算程序中實現起來相對比較容易，對于1階動力響應精確度較高。基于此，黏性邊界單元被嵌入到了很多大型商用軟件中以便于工程的實際應用，如LSDYNA、 FIAC 、ABAQUS[4]等。它的缺點是在施加人工邊界過程中只著重考慮了人工邊界輻射阻尼效應，忽略了人工邊界處半無限空間介質的彈性恢復性能。因此在實際工程應用中容易發生整體漂移的現象，進而使計算結果精度不高。透射邊界的優點是對于2階動力響應計算結果存在較高的精度，其缺點是動力計算分析程序的代碼編寫較復雜、在實際工程應用中可能存在高頻失穩的現象。黏彈性邊界不僅能有效的模擬半無限空間介質的輻射阻尼效應，同時還能較好的模擬半無限地基的彈性恢復性能。此外還具有良好的頻率穩定性和較高的精度。因此在實際工程中得到了廣泛的應用[5-10]。

作為大型通用有限元分析軟件之一的ANSYS軟件，雖然其自身已經存在著豐富的單元庫，但其單元庫中仍然存在著動力學求解中所需的單元，例如：黏彈元邊界單元、非線性單元、損傷單元，為了解決這一關鍵性問題，本文基于ANSYS平臺，采用FORTRAN語言開發了二維黏彈性邊界單元用戶子程序VSB_UELf，并結合ANSYS二次開發工具UPFs的功能特性，將其成功地嵌入到ANSYS中。最后，通過典型算例驗證了用戶自定義二維黏彈性邊界單元的正確性與可靠性，其計算結果的精度及穩定性均令人滿意。

1 二維黏彈性邊界單元

為了提高二維黏彈性邊界在土-結構相互作用求解過程中的計算精度，對于剛度系數KB和阻尼系數CB的選取顯得尤為重要，為了確定合理的KB和CB，其通常需要與人工邊界處的材料參數建立起相應的表達式[11-12]，其表達式如下：

(1)

切向阻尼系數：CBt=ρcs

(2)

(3)

(4)

為了便于對表達式中相關參數的理解，現對上式中的相關參數做詳細解釋，黏彈性人工邊界中的切向參數由αt表示，黏彈性人工邊界中的法向參數由αn表示。為了保持黏彈性邊界在地震動求解過程中的穩定性，依據αn和αt的推薦取值范圍，文中αt、αn的取值分別為0.5,1.0。其中二維黏彈性中的αn、αt的推薦取值詳見參考文獻[13]。

二維黏彈性邊界單元的剛度矩陣為：

(5)

二維黏彈性邊界單元的阻尼矩陣為：

(6)

2 用戶自定義二維黏彈性邊界單元的開發流程

在ANSYS軟件中的二次開發工具主要有4個：APDL、UPFs、UIDL、TclTk，其中UPFs是ANSYS二次開發的核心工具，UPFs是基于FORTRAN程序的二次開發，其主要功能是實現對ANSYS功能的擴充，即對ANSYS單元庫的擴充。同時ANSYS軟件也為用戶提供了諸多子程序，其中，UserElem.f子程序不僅可以幫助用戶傳遞自定義單元所需的的必要數據，同時還可以在動力計算過程中自動更新標準ANSYS中的數據庫和相關文件。需要說明的是UserElem.f子程序是基于ANSYS單元水平上的接口。基于UserElem.f子程序，用戶可以幾乎創建各種單元類型而無需直接訪問ANSYS數據庫和文件。因此,基于二維黏彈性人工邊界有關理論,借助于ANSYS二次開發工具中的核心工具UPFs以及接口子程序UserElem的功能特性，用戶自定義二維黏彈性人工邊界的開發具體流程圖見圖1。

為了實現將用戶自定義的二維黏彈性人工邊界單元成功的嵌入到ANSYS中，方便于推廣與應用，其中的關鍵環節之一是需要將接口子程序UserElem與標準程序進行成功連接，基于此用戶自定義的開發單元即可完成，其具體連接過程詳見參考文獻[14]。考慮到用戶自定義單元無法參與有限元單元的網格劃分，故應先通過ANSYS單元庫中的MESH200單元對所建立的二維有限元分析模型進行有限元網格剖分，最后再借助于ANSYS主程序中的內置命令將MESH200單元替換為對應的用戶自定義二維黏彈性單元即可。

圖1 二維黏彈性邊界單元開發流程圖

3 地震波的輸入方法

在對結構-地基動力相互作用問題進行分析時，為了在黏彈性人工邊界處實現地震波的輸入，借助于劉晶波等[15]介紹的一種地震波的直接輸入方法，即在邊界節點處采用力學中脫離體的概念，其主要思想是將輸入地震動轉化為作用于人工邊界節點上的等效荷載的方法來實現地震波動輸入。人工邊界節點處的等效荷載可表達為：

(7)

4 算例驗證

4.1 均質半空間二維自由場算例

如圖2所示均質二維彈性半空間，輸入荷載方程見式(8)，荷載按照第3節介紹的方法輸入，從均質二維彈性半空間的底部垂直輸入，進而獲得二維均質彈性半空間計算模型中典型位置的位移響應。

(8)

借助于大型通用軟件ANSYS平臺，在ANSYS中建立二維均質彈性半空間模型，該二維均質彈性有限域模型的計算尺寸：二維平面在長度方向的尺寸為40 m，寬度方向的尺寸為20 m。計算模型的網格剖分采用ANSYS單元庫中的四邊形單元Mesh200進行有限元網格劃分，剖分后的網格尺寸為每份1.0 m，剖分后的網格單元數及節點數分別為880 864。計算模型的相關材料參數如表1所示,網格剖分后的有限元計算模型如圖2所示，荷載輸入總持時為2 s，每一步的荷載時間步長取為0.005 s。

表1 二維均質彈性半空間計算模型材料參數

圖2 均質半空間二維計算模型網格剖分示意圖

為了考察在ANSYS中用戶自定義單元-二維黏彈性邊界單元的準確性與可靠性，選取計算模型中典型位置作為觀測點，五個觀測點的坐標分別為點A(-20，0)，B(0，0)，C(20，0)，D(20，-20)和E(0，-20)，圖3為各個觀測點在動荷載作用下的位移時程響應。從觀測點的位移時程曲線可以看出各個觀測點的位移時程曲線變化趨勢基本一致，最大偏差為1.5%；將A、B、C三點與D、E兩點的位移時程響應曲線對比來看，可以看出A、B、C三個觀測點的位移響應峰值接近D、E兩觀測點位移峰值的2倍，說明位移時程響應是由兩部分作用共同引起，即入射波與反射波共同作用。這些結論均與理論解基本相吻合，因此可以說明在ANSYS中用戶自定義單元——二維黏彈性邊界單元的可靠性與準確性。

圖3 A點—E點x向位移時程

4.2 層狀均質半空間二維自由場算例

為了考察在均質層狀土條件下的準確性與可靠性，借助于大型通用軟件ANSYS平臺，在ANSYS中建立了二維層狀均質彈性半空間模型(見圖4)，二維層狀均質彈性半空間計算模型在長邊方向上與寬度方向上的尺寸與與4.1節算例中模型的計算尺寸相同，并將計算模型分為兩層來進行計算，兩土層在ANSYS中用戶自定義單元——二維黏彈性單元邊界處的材料參數如表2所示,網格剖分后的有限元計算模型如圖4所示。

表2 二維層狀均質彈性半空間模型材料參數

(9)

式中：G(τ)=τ3H(τ)，τ=t/T，H(τ)為Heaviside函數分別按照黏彈性邊界及黏性邊界對模型進行計算，并與解析解[16]對比，觀測點的豎向位移響應對比結果如圖5所示。

從兩觀測點的豎向位移響應曲線(見圖5)可以看出，觀測點A與觀測點B的豎向位移響應曲線變化趨勢基本一致，僅位移峰值有一定的幅度差，A點的位移峰值較B點位移峰值所有減小，減小幅度為58%，其原因在于較波源的距離而言，B點比A點距離波源更近。此外，將三種不同解對比分析來看，黏彈性邊界解的變化趨勢及峰值與解析解的變化趨勢及峰值基本一致，說明黏彈性邊界計算結果具有高度的準確性與可靠性，而黏性邊界解較解析解以及黏彈性邊界解而言，在0.5 s前其變化趨勢與峰值與其它兩種解基本一致，而在0.5 s后出現了明顯的整體漂移現象，進而導致計算結果失真。如上所述，黏彈性邊界在考慮了無限半空間介質的彈性恢復性能后，很好地模擬了地震波在人工邊界處能量的傳播過程，同時也證明了在ANSYS中用戶自定義單元——二維黏彈性邊界單元具有一定的可靠性與準確性。

圖4 層狀半空間二維計算模型網格剖分示意圖

圖5 觀測點處豎向位移響應時程曲線

4.3 地震動輸入驗證

為了提高地震波以等效荷載的輸入形式計算結果可對比性，將考慮行波效應后的入射波位移時程作為該計算模型地表位移時程的理論解(見圖6)。從觀測點B的水平位移時程響應(見圖7、圖8)可以看出在兩種不同波即P波、S波的入射下，觀測點B的水平位移時程響應曲線的變化趨勢及位移峰值基本一致。同時得到的觀測點B水平位移時程響應與理論解的位移時程響應基本完全吻合，進一步可以說明，在ANSYS中用戶自定義單元——二維黏彈性邊界單元實現了在單元級別上的地震動等效荷載的施加，其施加方法及計算結果還具有較高的可靠性與準確性。

圖6 入射波位移時程

圖7 S波入射下B點水平位移響應

圖8 P波入射下B點水平位移響應

5 結 論

(1) 通過4.1節、4.2節的計算結果可以看出，借助于大型通用有限元有限元分析軟件-ANSYS-平臺，基于Fortran語言開發的用戶單元子程序VSB_UELf，在均質地基條件下及層狀土地基條件下均表現出計算結果的穩定性，同時計算精度也能滿足要求，說明了本文所開發的二維黏彈性邊界單元具有較好的準確性與可靠性。

(2) 通過4.3節的計算結果可以看出，在P波及S波的輸入下，計算結果均有著較好的穩定性與準確性，進而說明了本文所開發的二維黏彈性邊界單元中波動輸入方法是準確與可靠的。

(3) 另外，本文所開發的二維黏彈性邊界單元很容易擴展至三維，對于實際工程中的土-結構相互作用問題的求解有一定的使用價值，能更好的方便于實際工程的應用。