尹訓(xùn)強,羅 勇,王桂萱

( 大連大學(xué) 土木工程技術(shù)研究與開發(fā)中心，遼寧 大連 116622 )

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

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

1 二維黏彈性邊界單元

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

(1)

切向阻尼系數(shù)：CBt=ρcs

(2)

(3)

(4)

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

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

(5)

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

(6)

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

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

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

圖1 二維黏彈性邊界單元開發(fā)流程圖

3 地震波的輸入方法

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

(7)

4 算例驗證

4.1 均質(zhì)半空間二維自由場算例

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

(8)

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

表1 二維均質(zhì)彈性半空間計算模型材料參數(shù)

圖2 均質(zhì)半空間二維計算模型網(wǎng)格剖分示意圖

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

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

4.2 層狀均質(zhì)半空間二維自由場算例

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

表2 二維層狀均質(zhì)彈性半空間模型材料參數(shù)

(9)

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

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

圖4 層狀半空間二維計算模型網(wǎng)格剖分示意圖

圖5 觀測點處豎向位移響應(yīng)時程曲線

4.3 地震動輸入驗證

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

圖6 入射波位移時程

圖7 S波入射下B點水平位移響應(yīng)

圖8 P波入射下B點水平位移響應(yīng)

5 結(jié) 論

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

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

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