苑舉衛(wèi) 杜成斌 陳燈紅

（河海大學工程力學系,南京 210098）

在結(jié)構(gòu)-地基動力相互作用分析中為模擬近場能量向無限域的散射,通常需要在近場有限離散模型的人工截斷處設(shè)置人工邊界,以正確反映結(jié)構(gòu)-地基系統(tǒng)的動力特性.其中局部人工邊界由于其時空解耦特性而受到廣泛重視,目前較常用的有粘性邊界[1]、粘彈性邊界[2-4],以及廖振鵬[5]提出的透射邊界等.粘性邊界只考慮了邊界的耗能作用,而未考慮介質(zhì)的彈性恢復(fù)作用,在應(yīng)用中容易發(fā)生結(jié)構(gòu)整體飄移,精度較低；透射邊界雖具有較高精度,但在實際應(yīng)用中一般僅限于二階精度以內(nèi),并且編程較復(fù)雜、計算中可能引起高頻失穩(wěn)；粘彈性邊界能同時模擬波的散射效應(yīng)和半無限地基的彈性恢復(fù)能力,且能克服粘性邊界引起的低頻漂移問題,穩(wěn)定性好.

在考慮地基輻射阻尼的情況下,將地震動引入計算模型存在振動輸入和波動輸入兩種方式,前者適用于結(jié)構(gòu)剛度相對較小、輸入波頻率低、結(jié)構(gòu)尺寸遠較其需考慮的最短波長為小的情況,目前一般工業(yè)和民用建筑物和構(gòu)筑物多采用這種方式；而對于重要的大壩工程,由于壩體尺寸和重量都很大,地基對結(jié)構(gòu)體系動態(tài)特性的影響較大,這種情況一般都作為開放系統(tǒng)的波動問題求解,即將地震動輸入作為外源問題引入計算區(qū)域,波場分解法和等效邊界力法基于不同的思想均有效的建立了適用于粘彈性邊界的外源輸入方法[6-7].

ABAQUS作為功能最強大的有限元軟件之一,可以分析各種固體力學、結(jié)構(gòu)力學系統(tǒng),尤其是能夠處理非常復(fù)雜的問題和模擬高度非線性問題.ABAQUS還為用戶提供了方便靈活的二次開發(fā)平臺,包括若干用戶子程序（user subroutines）以及實用程序（utility routines）,為解決實際工程的復(fù)雜問題提供了方便.

為方便地實現(xiàn)粘彈性人工邊界的設(shè)置及地震動的有效波動輸入,本文利用有限元軟件ABAQUS提供的自定義單元子程序UEL結(jié)合隱式積分算法[8],通過定義邊界單元對整個系統(tǒng)的雅克比矩陣貢獻來實現(xiàn)粘彈性人工邊界.基于波場分離方法,將人工邊界上的總波場分解為無局部場地效應(yīng)影響的自由場與局部場地效應(yīng)引起的散射場,散射場由粘彈性人工邊界吸收,則作用于人工邊界的無局部場地效應(yīng)影響的自由場可由設(shè)計地震動位移時程和速度時程統(tǒng)一表達,利用子程序DSLOAD、UTRACLOAD描述邊界的加載歷程.將子程序 UEL、DSLOAD、UT RACLOAD以及設(shè)計地震動等數(shù)據(jù)匯集成 SDAB.FOR,結(jié)合在ABAQUS中建立的有限元模型,可方便有效地進行結(jié)構(gòu)-地基相互作用的分析.

1 粘彈性人工邊界的有限元實現(xiàn)

粘彈性人工邊界描述的是由近場向遠域散射的外行波在人工截斷邊界處滿足的應(yīng)力條件,可表示如下

式中,σ表示邊界點的應(yīng)力向量；u、˙u分別為邊界點的位移和速度向量；KB、CB分別為以粘彈性人工邊界彈簧系數(shù)、粘性系數(shù)為元素的對角矩陣.粘彈性人工邊界的實質(zhì)是在人工邊界上施加切向和法向方向上的彈簧和阻尼器,其系數(shù)按式（2）計算

式中,KBN、KBT分別表示法向和切向彈簧系數(shù)；CBN、CBT分別為表示法向和切向阻尼系數(shù)；rb為散射波源至人工邊界的距離；cp、cs分別為P波和S波波速；G為介質(zhì)剪切模量；λ為拉梅常數(shù)；ρ為介質(zhì)質(zhì)量密度；α、β為無量綱參數(shù),分別取0.8、1.1.

粘彈性人工邊界作為一種應(yīng)力連續(xù)分布的人工邊界條件,可采用與普通單元類似的形函數(shù)進行離散,即一致粘彈性人工邊界,經(jīng)離散后,粘彈性人工邊界單元的剛度矩陣[9]如式（3）,阻尼矩陣具有相同的形式,只需將剛度系數(shù)K用相應(yīng)的阻尼系數(shù)C替換即可.

在計算區(qū)域邊界外側(cè)拓展一層六面體單元,單元外側(cè)節(jié)點施加固定約束,一致粘彈性人工邊界單元的剛度矩陣、阻尼矩陣通過ABAQUS的自定義單元子程序UEL,結(jié)合動力分析所采用的 Hilber-Hughes-T aylor隱式積分方法,以邊界單元對整個系統(tǒng)的矩陣貢獻的形式寫入.單元的厚度只有象征意義,沒有限制,對結(jié)果沒有影響.

其中,Sb為三維一致粘彈性人工邊界的面積[10].

2 地震動輸入方式

2.1 粘彈性邊界外源波動輸入

基于波場分離,將人工邊界上的總波場分解為無局部場地效應(yīng)影響的自由場與局部場地效應(yīng)引起的散射場,局部場地效應(yīng)引起的散射場由粘彈性人工邊界條件可自動滿足,將外源地震波引入計算區(qū)域,只需在人工邊界處施加自由場荷載,把近域地基作為半無限空間自由場地基的子結(jié)構(gòu),由其在入射地震波作用下的邊界相互作用力[11]給出:

式中,n為邊界外法線方向余弦向量；KB、CB均與式（1）含義相同分別為自由場位移矢量 、速度矢量和應(yīng)力矢量,當在近域地基底邊邊界輸入的地震波波陣面為平面時,自由場應(yīng)力矢量可由速度場表示.

假定地震波為垂直入射的平面P波和平面S波,則根據(jù)波動傳播規(guī)律及波場應(yīng)力狀態(tài)邊界相互作用力（4）可由設(shè)計地震動位移時程和速度時程統(tǒng)一表達,以分布荷載的形式施加于人工邊界,其中邊界法向荷載采用子程序DSLOAD描述,切向荷載由子程序UTRACLOAD描述.

2.2 粘彈性邊界內(nèi)源地震動輸入

在計算區(qū)域人工截斷處設(shè)置粘彈性人工邊界來反映無限區(qū)域的情況下,地震荷載的引入還有另一種輸入方式,即考慮粘彈性邊界的內(nèi)源地震動輸入方法.在均勻輸入情況下,或者空間不均勻地震動較小的情況下,地震動產(chǎn)生的荷載只與結(jié)構(gòu)質(zhì)量有關(guān),而與地基質(zhì)量無關(guān),但是系統(tǒng)的響應(yīng)包含地基質(zhì)量,這也是與常規(guī)無質(zhì)量地基不同的地方[12].

因此,本文將考慮3種地震動輸入方式,即考慮粘彈性邊界條件的外源波動輸入（輸入方式1）和內(nèi)源輸入（輸入方式2）,及常規(guī)的無質(zhì)量地基輸入（輸入方式3）,比較不同輸入方式對結(jié)構(gòu)-地基系統(tǒng)動力響應(yīng)的影響.

3 外源波動算例驗證

考察均勻三維彈性半空間,從底部人工邊界分別豎直入射平面S波、P波的動力響應(yīng).入射位移脈沖波方程為

從均勻三維彈性半空間截取長和寬均為400m、高度為600m的有限區(qū)域作為計算區(qū)域,用有限單元法進行空間離散,3個方向網(wǎng)格尺寸均為20 m,共離散為13671個結(jié)點,12000個六面體單元.在計算區(qū)域外側(cè)及底側(cè)拓展一層六面體單元,構(gòu)造一致粘彈性邊界單元.自由場荷載作用于人工截斷處各表面.材料參數(shù)取為:彈性模量E=4.88GPa,密度ρ=2000 kg/m3,泊松比ν=0.22.P波波速cp=1669.05m/s,S波波速cs=1000m/s,時間步長取為Δt=0.01 s.

該問題的自由地表位移解析解為考慮行波延遲后放大2倍的入射位移時程.圖1、圖2分別為各觀測點在平面S波入射時的水平向位移響應(yīng)和P波入射時的豎向位移響應(yīng).可見入射波到達地表后與反射波疊加,地表的響應(yīng)為入射波的2倍；底面觀測點在經(jīng)歷入射波和經(jīng)地表反射的波后響應(yīng)為 0,波傳向遠域,均與解析解一致.數(shù)值模擬結(jié)果符合波的傳播規(guī)律,且有良好的模擬精度,說明本文的三維粘彈性邊界單元及外源波動輸入子程序SDAB.FOR是準確的.

圖1 S波入射時底面點和地表點的位移時程

圖2 P波入射時底面點和地表點的位移時程

4 工程應(yīng)用

某水電站位于我國瀾滄江上游河段,壩體為拋物線變厚度雙曲拱壩,壩底高程953 m,壩頂高程1245 m,最大壩高292m,壩頂弧長 935m,壩頂厚 12m,壩底厚73m,壩體混凝土動態(tài)彈性模量Ed=20GPa,密度 ρd=2400kg/m3,泊松比 νd=0.189,地基彈性模量Ed=27.3 GPa,密度 ρd=2400 kg/m3,泊松比 νd=0.25.

根據(jù)場地地震安全性評價成果,工程區(qū)設(shè)計烈度為Ⅸ度,水平向設(shè)計地震動峰值加速度為0.308g,豎向設(shè)計地震動峰值加速度取水平向的2/3,即0.205 g,根據(jù)規(guī)范譜人工合成地震波如圖3所示,計算總時間為10 s,時間步長為0.005s.

采用文獻[13]對不同地基范圍粘彈性人工邊界的精度研究結(jié)論,地基部分沿結(jié)構(gòu)各方向延伸一倍壩高就可滿足工程精度要求,邊界設(shè)置粘彈性人工邊界.壩體和基巖用三維8結(jié)點六面體單元進行有限元離散,共劃分31026個單元,35366個結(jié)點,有限元網(wǎng)格如圖4所示.圖5給出了個關(guān)鍵部位的位置示意圖,括號中表示對應(yīng)下游面.

圖4 拱壩-地基系統(tǒng)有限元模型 圖5 關(guān)鍵部位示意圖

在上述3種不同地震動輸入方式下,上、下游面關(guān)鍵部位的應(yīng)力峰值如表1、2所示,其中括號中百分數(shù)表示與常規(guī)地震動輸入方法計算結(jié)果的降低幅度,拱冠處上下游的主拉應(yīng)力時程曲線如圖6所示.從表中可以看出,除個別點外,采用無質(zhì)量地基模型的地震動常規(guī)輸入方法的響應(yīng)最大,引入粘彈性人工邊界后,壩體的動力響應(yīng)顯著降低,且采用外源波動輸入方法計算的響應(yīng)最小.

表1 上游面各關(guān)鍵部位應(yīng)力比較 （單位:MPa）

續(xù)表1 上游面各關(guān)鍵部位應(yīng)力比較（單位:MPa）

表2 下游面各關(guān)鍵部位應(yīng)力比較 （單位:MPa）

外源波動輸入方法與常規(guī)地震動輸入方法相比,無論是主拉應(yīng)力還是主壓應(yīng)力降幅都較大,最大降幅達53%.相互作用無質(zhì)量地基與常規(guī)地震動輸入方法相比,主拉應(yīng)力降幅較大,其中上游面主拉應(yīng)力降幅在10%～35%之間,下游面主拉應(yīng)力降幅在15%～45%之間,上下游面主壓應(yīng)力降幅均25%.這兩種輸入方法相比在拱冠處的響應(yīng)較其他部位接近.

5 結(jié) 論

本文利用FORT RAN語言編寫了適用于粘彈性人工邊界的外源波動輸入子程序SDAB.FOR,包含粘彈性邊界和加載兩部分,從而可方便地實現(xiàn)地震動的有效輸入.將考慮粘彈性邊界條件的外源波動輸入和內(nèi)源地震動輸入方法分別應(yīng)用于某高拱壩的動力響應(yīng)分析,并與常規(guī)無質(zhì)量地基模型計算結(jié)果進行了對比,結(jié)果表明,考慮地基輻射阻尼的兩種地震動輸入方法均會顯著的降低結(jié)構(gòu)的動力響應(yīng),其中采用外源波動輸入方法的降幅更大.

[1]Lysmer J,Kulemeyer R L.Finite Dynamic Model for Infinite Media[J].Journal of Engineering Mechanics Division,ASCE,1969,95（4）:759-877.

[2]Deeks A J,Randolph M F.Axisymmetric Time-domain Transmitting Boundaries[J].Journal of Engineering Mechanics Division,ASCE,1994,120（1）:25-42.

[3]劉晶波,谷 音,杜義欣.一致粘彈性人工邊界及粘彈性邊界單元[J].巖土工程學報,2006,28（9）:1070-1075.

[4]杜修力,趙 密,王進廷.近場波動模擬的一種應(yīng)力人工邊界[J].力學學報,2006,38（1）:49-56.

[5]廖振鵬.工程波動理論導(dǎo)論[M].北京:科學出版社,2002.

[6]劉晶波,呂彥東.結(jié)構(gòu)-地基動力相互作用問題分析的一種直接方法[J].土木工程學報,1998,31（3）:55-64.

[7]趙建鋒,杜修力,韓 強等.外源波動問題數(shù)值模擬的一種實現(xiàn)方式[J].工程力學,2007,24（4）:52-58.

[8]Abaqus Theory Manual（Version 6.7）,Dassault Systems Inc,2007.

[9]谷 音,劉晶波,杜義欣.三維一致粘彈性人工邊界及等效粘彈性邊界單元[J].工程力學,2007,24（12）:31-37.

[10]徐 磊,葉志才,任青文等.基于ABAQUS的粘彈性動力人工邊界精確自動施加[J].三峽大學學報:自然科學版,2010,32（1）:20-23.

[11]陳厚群.壩址地震動輸入機制探討[J].水利學報,2006,37（12）:1417-1423.

[12]李 靜,陳健云.動力相互作用分析中無質(zhì)量地基的應(yīng)用研究[J].世界地震工程,2007,23（2）:58-62.

[13]Zhang Chuhan,Pan Jianwen,Wang Jinting.Influence of seismic input mechanisms and radiation damping on arch dam response[J].Soil Dynamics and Earthquake Engineering,2009,29（9）:1282-1293.