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

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

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

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

1 二維黏彈性邊界單元

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

(1)

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

(2)

(3)

(4)

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

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

(5)

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

(6)

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

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

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

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

3 地震波的輸入方法

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

(7)

4 算例驗(yàn)證

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

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

(8)

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

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

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

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

圖3 A點(diǎn)—E點(diǎn)x向位移時(shí)程

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

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

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

(9)

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

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

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

圖5 觀測(cè)點(diǎn)處豎向位移響應(yīng)時(shí)程曲線

4.3 地震動(dòng)輸入驗(yàn)證

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

圖6 入射波位移時(shí)程

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

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

5 結(jié) 論

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

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

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