多層波導中矢量波動的時域人工邊界條件1)

吳利華 趙 密 杜修力

(北京工業大學城市與工程安全減災教育部重點實驗室,北京 100124)

引言

用有限元法等數值方法分析大型結構的瞬態動力響應,一般的做法是引入人工邊界,將結構?無限地基系統劃分為有限域和無限域.無限域被截去,有限域用有限元法等數值方法模擬.為了消除波在人工邊界上的虛假反射,需要在有限域的人工邊界上施加人工邊界條件(artificial boundary condition,ABC)來描述波動在無限域中的傳播[1].ABC 是基于無限域的控制方程推導得到的.一般假定無限域是線性的,因為通常認為無限域中只存在向無窮遠傳播的波,不存在向有限域傳播的波.

剛性基巖上臥成層土是一種普遍存在的地基形式,又稱其為半開放的波導模型.與開放的半空間模型中只存在傳向無窮遠的行波相比,波導模型由于截止頻率的存在,介質中同時存在截止頻率以下的近場快衰波和截止頻率以上的傳向無窮遠的行波[2-3].對于半空間模型,行波的振幅隨著傳播距離的增大快速減弱,所以即使不考慮介質材料阻尼的吸收效應,其無限域也可以被假定為小應變,即可以被假定為線性的.而對于波導模型,由于能量被困在介質中,行波的振幅衰減得較為緩慢,如果僅考慮幾何散射不考慮介質吸收效應,將無限域假定為小應變是有問題的,所以更合理的做法應該是將波導模型考慮為有材料阻尼的線彈性介質.

本文的研究對象是具有確定尺寸參數的多層波導結構.Wang 等[4]研究表明,當多層波導內部結構尺寸存在隨機缺陷時,彈性波會產生局部化現象.該現象使得在原本可以傳播的頻域范圍內彈性波被禁止傳播,進而多層波導結構中會產生不均勻的應力分布.鑒于問題的復雜性以及篇幅的限制,本文的研究范圍不包括上述多層波導中隨機失諧引起的彈性波局部化問題.

如有限域含有非線性等許多情況下,需要ABC是時域方法.相較于半空間模型的時域ABC 已有大量的研究工作[5-13],波導模型的時域ABC 研究工作較少.由于波導模型的ABC 需要能夠同時模擬快衰波和行波的能量傳遞,直接在時域下建立其ABC 是有挑戰的.根據是否精確地處理波導模型無限域定解問題,以下將已有的波導模型的ABC 分為精確的ABC 和近似的ABC 進行闡述.

精確的ABC 一般是首先基于解析法或半解析半離散法獲得精確的頻域動力剛度,再通過時域化方法建立時域ABC.解析法一般對單層介質標量波[14-16]適用,多層介質的解析解難以獲得.半解析半離散法只對無限域的人工邊界進行離散,對于單層介質和多層介質都適用,如邊界元法(BEM)[17]、薄層法(TLM)[18-21]和比例邊界有限元法(SBFEM)[3,22-26].BEM 需要基本解,而多層介質格林函數的基本解是復雜的,基本解的存在與否限制了BEM 在工程實際中的應用.相比于BEM,TLM 和SBFEM 不需要基本解,能夠自動滿足無窮遠輻射邊界條件.將頻域動力剛度變換到時域,一種方法是時間卷積計算[3,25],該方法不僅耗時,而且不適用于非線性分析.學者們通過時間局部化方法建立了精確穩定的時間局部的ABC.Zhao 等[15]基于解析法獲得2D 單層波導標量波的動力剛度,再用有理函數逼近方法建立精確的時域ABC.Prempramote 等[22]將解析法獲得的動力剛度用雙漸近連分式展開,建立了2D 單層標量波的時域ABC.Li 等[16]基于解析法獲得2D 和3D 均勻層標量波的動力剛度,用動力剛度連分式展開的時間局部化方法建立了精確的時域ABC.Liu 等[26]基于半解析半離散法,利用算子分離法獲得2D 單層標量波的時域ABC.Wang 等[23]基于SBFEM 和動力剛度的雙漸近連分式展開[22]建立了2D 庫水的時域ABC,并將方法嵌入到通用有限元軟件ABAQUS 中用來計算壩?庫水的動力相互作用.但是上述精確的時域ABC 都是針對單層波導標量波獲得.根據作者的調查,目前幾乎沒有關于多層波導穩定的精確時域ABC 的相關報導.因為對于多層波導,其精確的時域ABC 易發生失穩,目前還沒有較好地抑制其失穩的方法.

多層波導精確的時域ABC 難以獲得,學者們通過對無限域定解問題作適量的簡化,建立了多層波導近似的時域ABC.近似的低階時空局部的ABC 由于其形式簡單,在工程中被廣泛使用,如黏性邊界[5,27]、多次透射邊界[7,11]以及黏彈性邊界[9-10].這類低階時空局部ABC 是基于半空間模型建立的,用于波導模型時精度較低.Hagstrom 等[28]建立了多層波導標量波的高階時空局部ABC.完美匹配層法[29-30]是近年的研究熱點,它是通過在人工邊界外附加一層具有耗能作用的緩沖區來吸收透射波,但吸收層具有一定厚度,層厚度和耗能參數選擇不當時,會導致精度低或發生失穩.Zhao 等[32]基于半解析半離散法和一種全頻域收斂的動力剛度連分式展開構建了多層波導標量波的時域ABC,并對其穩定性進行了證明.吳利華等[33]將多層波導考慮為達朗貝爾黏彈性介質,利用和文獻[32]類似的思路建立了考慮阻尼的多層波導標量波的時域ABC.高毅超等[34]將彈性動力學矢量方程強行解耦為兩個標量波方程,基于SBFEM和雙漸近連分式[22]構建了一種模擬單層波導矢量波的近似ABC.Liu 等[35]將文獻[34]的方法嵌入到開源有限元軟件OpenSees 中,分析了地下車站的地震動響應.Liu 等[31]基于半解析半離散法和算子分離法建立了多層波導矢量波的近似ABC,但該方法有時會出現失穩現象.就作者所知,目前幾乎沒有專門針對多層波導矢量波的穩定時域ABC 的相關研究.

本文的目標是建立多層波導矢量波穩定的時域ABC.借鑒文獻[34]的做法將多層波導矢量波動方程強行解耦,得到x方向和y方向的兩個標量波動方程.將文獻[33]構建達朗貝爾黏彈性多層波導標量波的時域ABC 的思路用于解耦的x方向和y方向的兩個標量波動方程,構建含有瑞利阻尼的線彈性多層波導矢量波的穩定的時域ABC.

1 問題描述

用有限元法分析如圖1 所示的剛性基巖上臥成層地基上結構在動力載荷作用下平面內的時程響應.由于有限元法只能分析有限區域,需要引入人工邊界將左右兩個半無限域截去.為了避免波在人工邊界處引起虛假反射,需要在有限域的人工邊界處施加一個人工邊界條件(ABC)來模擬波在無限域中的傳播.

圖1 剛性基巖上臥成層土Fig.1 A multilayered medium overlying rigid bedrock

有限域可以是非線性的,其有限元方程為

其中,下標b 和i 分別表示人工邊界部分和有限域除了人工邊界的其他部分,M是質量矩陣,C是阻尼矩陣,K是剛度矩陣,u是平面內位移向量,是非線性內力,b是體力,fi是施加在有限域的載荷,fb是施加在有限域人工邊界上的力,即本文要建立的ABC,用來表示被截掉的無限域和有限域的相互作用.

一般認為無限域是線性的,若無限域包含非線性,需要對其進行等效線性化處理,使其滿足線性條件.無限域的波動方程為

其中,變量上方的點表示對時間求導,(x,y)表示笛卡爾坐標,上標T 表示轉置;U={ux uy}T,ux,uy分別為x方向和y方向的位移,ux和uy在不同土層的分界面以及人工邊界處是連續的;ρ 是質量密度,λ 是拉梅常數,G是剪切模量,每一層的材料參數是常量,不同層的材料參數可不同.另外,無限域的上表面物理邊界條件是應力為0,下表面物理邊界條件是位移為0.土層的初始狀態為靜止.

2 時域人工邊界條件

將無限域作為研究對象,建立時域吸收邊界條件.文獻[33]建立了模擬剛性基巖上臥多層介質中標量波傳播的ABC,本文將該方法擴展到矢量波.方法的思路是:(1)將矢量波動方程簡化為x方向和y方向解耦的兩個標量波動方程;(2)基于比例邊界有限元法,得到兩個標量波動方程模態空間下半離散的頻域動力剛度方程;(3)通過文獻[33]中提出的連分式來表示頻域動力剛度;(4)引入輔助變量,將連分式時域化,得到時域吸收邊界條件;(5)所建立的ABC 可以和有限域的有限元方程無縫耦合,通過數值積分算法求解.

2.1 彈性波動方程簡化

高毅超等[34]通過將矢量波動方程強行解耦,基于雙漸近連分式建立了單層波導模型的ABC.本文借鑒其做法,忽略無限域波動方程(2)中x和y的耦合項,簡化后得到關于x方向和y方向的兩個標量方程分別為

其中,x和y解耦簡化處理后的應力?應變關系為

經過解耦簡化處理后的上下表面物理邊界條件分別為

其中,H為半無限域的總層高.

2.2 模態空間下的頻域動力剛度

根據比例邊界有限元法(SBFEM)[23]建立模態空間下的頻域動力剛度方程.

沿著土層y方向半離散簡化處理后的無限域波動方程(3),同時考慮簡化后的物理邊界條件式(4)及介質交界面的連續條件,得到無限域簡化彈性波動方程的比例邊界有限元(SBFE)方程

基于簡化后的本構關系得到無限域的人工邊界上x和y方向的等效結點力分別為

引言的第二段闡述了波導模型應該考慮阻尼的吸收效應.由于瑞利阻尼是目前應用最廣泛的阻尼矩陣的數值模型,本文假定無限域土層含有瑞利阻尼.瑞利阻尼定義阻尼矩陣是質量陣和剛度陣的線性組合[36],即CRayleigh=a0M+a1K,其中,ξ 為陣型阻尼比,ωi,ωj為選擇用于確定阻尼系數a0和a1的頻率點.瑞利阻尼的引入體現在式(5a)和式(5b)的第四項,需要說明的是,這里無限域瑞利阻尼的質量陣和剛度陣都是通過半離散x和y解耦的標量控制方程得到的,注意其與二維有限元法中的瑞利阻尼相區別.

為了使問題簡單化,利用模態分解法[36],將物理空間下的SBFE 方程變換到模態空間.計算如下廣義特征值問題

其中,上標?1 表示對矩陣求逆.

2.3 動力剛度的連分式展開

其中,z=iω/ωn,β=a0/ωn+a1ωn,i 是虛數單位,ω 表示圓頻率.x和y方向分別有c=cp和c=cs.

由于x方向和y方向的SBFE 方程形式相同,以下僅給出x方向的推導過程,y方向的推導同x方向.以下推導思路類似文獻[33],本文僅給出與文獻[33]不同的地方以及相關重要公式,詳細推導過程可參見文獻[33].根據式(12)和式(13),無限域x方向模態空間的動力剛度可表示為如下的連分式

和文獻[33]一樣,連分式(14)也可用來近似地表示多層介質的動力剛度.對于多層介質,連分式中的通過求解如下黎卡提方程得到

2.4 基于連分式的時域吸收邊界條件

同文獻[33]一樣,通過引入輔助變量的方式將連分式時域化.將式(14)代入式(10a),同時引入輔助變量,得到

如今，微課的設計絕大多數都是對教材中知識點的教授，還滯留在對學生采取知識灌溉的方式上，然而微課的作用絕不止在于此。微課的時間很短，但是使用微課的最終目標是希望做到事半功倍，以小見大。所以，教師必須從簡潔的授課型微課入手，漸漸向啟發型、探究型的教學模式發展。想要做到這一點，教師在使用微課之前，就要做好充實的資料收集與課前準備，以便能夠科學有效地權衡探究型活動和師生互動使用的微課時間。

式(20)中,人工邊界上單個結點的水平自由度和豎向自由度是空間解耦的,但是人工邊界結點間的水平自由度和結點間的豎向自由度都是空間耦聯的.其中

2.5 人工邊界條件與有限元法耦合

本文建立的ABC 可與有限元法無縫耦合.將式(20)代入有限元方程式(1),得

式(21)可通過標準的時間積分算法如Newmark 常平均加速度法等求解.

由于在有限域的人工邊界上施加ABC 增加的額外自由度為

其中,nL和nR分別表示左無限域和右無限域參與計算的模態數.J表示連分式的階數,通過算例分析統計,一般可取J=3,J的取值會在第3 節詳細闡述.對于實際地震工程等問題,根據經驗一般取模態數n≤3.所以額外增加的自由度數一般不會超過36,與有限域的自由度數相比,額外增加的自由度數幾乎可以忽略不計.

本文建立的時域ABC 是對無限域矢量波動的近似模擬,為了實現高精度,需要將有限域的人工邊界放置在距離興趣域L處,L的取值將在第3 節通過數值算例詳細分析.

本文的ABC 是基于全頻域收斂穩定的連分式建立的,所以理論上ABC 應該是時域穩定的.方法的穩定性將在第3 節通過數值算例進一步說明.

3 方法的性能

本節通過數值算例說明方法的精度和穩定性.本文方法包含3 個可變參數:無限域的模態數n,連分式階數J和人工邊界遠離興趣域的距離L.文獻[33]表明,僅用能被載荷激起的無限域的模態數n參與計算,不但能減少計算量,同時還不影響精度.下文詳細分析J和L的取值要求.

分析如圖2 所示的3 種模型.Model 1 和Model 2 都是自由場地,Model 1 的總層高H為20 m,Model 2 的總層高H為40 m.比較Model 1 和Model 2 的結果可以考察L,J取值與H的關系.Model 3 內含一個尺寸為5×5 的地下方洞,其總層高H為20 m.比較Model 1 和Model 3 的結果可以考察有無地下結構及地下結構尺寸對L,J取值的影響.3 種模型都是在地表作用10 m 長的水平線載荷,線載荷為狄拉克脈沖,其時程及頻譜如圖3 所示.為了考察狄拉克脈沖寬度是否會影響本文方法的精度和收斂性,考慮三種狄拉克脈沖:T=0.4,0.2 和0.1,其相應的頻譜范圍分別為0～10 Hz,0～20 Hz 和0～40 Hz.圖2 中每個模型的興趣域(region of interest,ROI)為藍虛線框起來的部分,紅實線表示人工邊界,兩條紅實線框起來的部分是用于計算的有限域.L表示人工邊界到興趣域的距離.

圖2 數值算例模型Fig.2 Three models for numerical examples

圖3 狄拉克脈沖(T=0.4,0.2 和0.1)Fig.3 Dirac pulse(T=0.4,0.2 and 0.1)

為了分析土層材料參數對L,J取值的影響,將圖2 中3 種模型分別都考慮為均勻介質和分別都考慮為4 層介質.分別將單層Model 1,單層Model 2 和單層Model 3 記為Case-1,Case-2 和Case-3;分別將4層Model 1,4 層Model 2 和4 層Model 3 記為Case-4,Case-5 和Case-6.Case-1～Case-6 土層都被考慮為有瑞利阻尼的線彈性介質,其土層材料參數見表1.6 種工況的計算時長都為2 s,輸出A 點的水平位移時程和水平加速度時程.

用有限元法分析Case-1～Case-6.圖2 所示的Model 1～Model 3 中兩條紅實線表示的人工邊界框起來的部分為有限域,人工邊界到興趣域的距離被標記為L,在人工邊界處施加ABC.有限域用四邊形網格離散,Newmark 常平均加速度法用于時間積分計算.狄拉克脈沖中T=0.4 時,有限域的網格尺寸為2.5 m×2.5 m,時間步長為0.005 s.T=0.2 和0.1 時的網格尺寸和時間步長分別是T=0.4 時的0.5 倍和0.25 倍.

表1 Case-1～Case-6 的土層材料參數及其瑞利阻尼系數Table 1 Material parameters and Rayleigh damping coefficients of Case-1～Case-6

用于確定瑞利阻尼系數的頻率點的選擇方式有較多研究[37].本節的重點是驗證本文提出方法的精度和穩定性,由于瑞利阻尼系數的具體數值不影響本文方法的性能,所以這里采用簡單的頻率點的選取方式,取土層前2 階頻率用于確定其瑞利阻尼系數a0和a1.Case-1～Case-6 土層的阻尼比ξ 都取為0.05,表1 列出了6 種工況的瑞利阻尼系數值.狄拉克脈沖中T=0.4 時,Case-1～Case-6 用于計算的無限域模態數n分別為2,3,2,1,3 和1.

將足夠大有限域的有限元結果作為驗證本文方法的參考解,大有限域尺寸的選取原則是人工邊界上的虛假反射波在計算時間內不會污染興趣域.

用數值算例結果分析J和L的取值對方法精度的影響.定義如下時程相對誤差

其中,Xj為基于有限元法和ABC 得到的小有限域的時程結果,Yj為僅基于有限元法得到的大有限域的時程結果(即參考解),N為總時步.

圖4 展示了Case-1～Case-6 在連分式階數J分別為1,2,3,4,5 的情況下,L=H～5H時對應的A點水平加速度時程的相對誤差,這里狄拉克脈沖中T=0.4.為了說明本文方法的優越性,圖4 也展示了同樣L下黏性邊界[5](VB)、劉晶波等[9]提出的黏彈性邊界(VSB-Liu)和杜修力等[10]提出的黏彈性邊界(VSB-Du)的結果.黏性邊界和黏彈性邊界都是經典的ABCs,因其形式簡單且時域穩定,在工程中被廣泛使用.黏性邊界[5]是基于平面波假定建立的阻尼器邊界條件,其每層的阻尼元件法向參數為CN=ρcp,切向參數為CT=ρcs.黏彈性邊界是基于柱面波假定建立的并聯彈簧?阻尼器邊界條件.劉晶波等[9]提出的黏彈性邊界每層的彈簧元件法向參數為KN=2G/r,切向參數為KT=3G/(2r),每層的阻尼元件法向參數為CN=ρcp,切向參數為CT=ρcs;杜修力等[10]提出的黏彈性邊界每層的彈簧元件法向參數為KN=(λ+2G)/(3.6r),切向參數為KT=G/(3.6r),每層的阻尼元件法向參數為CN=1.1ρcp,切向參數為CT=1.1ρcs.3 個公式中每層的材料參數按實際值取,r認為是有限域的豎向對稱軸到人工邊界的水平距離.從圖中可以看出,6 種工況呈現的結果規律類似.每種工況下,黏性邊界和2 種黏彈性邊界結果的相對誤差接近,都遠大于本文方法的結果.本文方法的J=1,2,3,4,5 對應的5 條曲線都是J=1 時相對誤差較大,J=2 時相對誤差快速變小,J=3,4 和5時三者的結果基本重合.根據大量算例結果(不僅限于圖4 所列的算例結果)統計,當J>3 時,再增加J的值,基本不再改善精度.這是因為本文方法是一種近似方法,連分式只是近似地表示無限域的動力剛度,即使取很高的階數,連分式也不能逼近無限域動力剛度的精確解.因近似損失的精度需要通過增大人工邊界到興趣域的距離L來彌補.為了減少可變因素,作者建議在使用本文方法時,參數J可以統一取為3.

圖4 6 種工況下基于本文方法、黏性邊界(VB)和2 種黏彈性邊界(VSB-Liu 和VSB-Du)的A 點水平加速度時程的相對誤差Fig.4 Relative errors of time-history horizontal acceleration solutions at point A from the proposed method,the viscous boundary and two kinds of viscous-spring boundary VSB-Liu and VSB-Du for Case-1～Case-6

進一步詳細分析L的選取原則.分別計算Case-1～Case-6,當J=3,L為H～5H時A點水平位移時程和水平加速度時程各自的相對誤差,這里狄拉克脈沖中T=0.4.將相對誤差≤5%時各自對應的L值,以及同樣L下黏性邊界和黏彈性邊界結果的相對誤差列于表2.根據表2 可以得出以下結論(表中第2 列u表示位移,a表示加速度;第3 列表頭L表示人工邊界到興趣域的距離;第4 列表頭ABC 表示本文方法的結果;第5 列表頭VB 表示黏性邊界[5]的結果;第6 列表頭VSB-Liu 表示劉晶波等[9]提出的黏彈性邊界的結果;第7 列表頭VSB-Du 表示杜修力等[10]提出的黏彈性邊界的結果).

(1)對比表2 的最后3 列可以看出,黏性邊界和黏彈性邊界結果的相對誤差接近,都約為本文方法結果的5～6 倍左右,除了Case-2 和Case-5 的加速度時程結果,其黏性邊界和黏彈性邊界結果的相對誤差是本文方法結果的2～3 倍.另外,本文方法與黏性邊界和黏彈性邊界相比,僅多出少量的輔助自由度數,對于Case-1～Case-6,其值分別為24,36,24,12,36,12.因此本文方法相比于黏性邊界和黏彈性邊界,能夠在幾乎不增加計算量的情況下使精度提高約2～6倍.精度的提高程度與土層的總厚度有關,原因是本文方法是針對層模型建立的,其應用在深厚土層和淺土層都具有高精度;黏性邊界和黏彈性邊界是針對半空間模型建立的,由于深厚土層相較于淺土層更接近半空間模型,因此其在深厚土層中要比在淺土層中的精度高.

(2)對比表2 第3 列表示的水平位移對應的L值和水平加速度對應的L值可以看出,每種工況下,水平加速度對應的L值都是水平位移對應結果的2 倍,除了Case-5 是2.5 倍.說明人工邊界的位置以加速度為控制指標來確定更為嚴格.

(3)比較單層的Case-1 和Case-3 水平加速度對應的L值.Case-1 和Case-3 總層高相同,Case-1 為自由場地,Case-3 含有5 m×5 m 的地下結構.兩者都是L=3H時,水平加速度時程的相對誤差能控制在5%以內.說明L的取值基本不受是否含有地下結構或結構尺寸影響.比較4 層的Case-4 和Case-6 能得出同樣的結論.

(4)比較單層Case-1 和Case-2 水平加速度對應的L值.Case-1 和Case-2 都是自由場地,Case-2 的總層高H是Case-1 的2 倍.兩者都是L=3H時,水平加速度時程的相對誤差能控制在5%以內.同樣,比較4 層的Case-4 和Case-5.為了將水平加速度時程的相對誤差控制在5%以內,Case-4 要求L=2H,Case-5要求L=2.5H.說明L的取值大約與土層總層高H成正比關系,關系系數與土層的材料參數有關.根據大量算例統計,一般L=3H都能將加速度時程相對誤差控制在5%以內.

表2 本文方法的A 點水平位移時程和水平加速度時程各自的相對誤差≤5%時對應的L 值,以及同樣L 下黏性邊界和2 種黏彈性邊界結果的相對誤差Table 2 Values of L when relative errors of horizontal displacement and of horizontal acceleration at point A from the proposed method not greater than 5%,and relative errors of those solutions from the viscous boundary and two kinds of viscous-spring boundary under the same value of L

圖5 展示了本文方法的計算結果滿足工程精度要求(相對誤差≤5%)時的A點水平加速度時程.其中,J=3,Case-1～Case-6 對應的L值分別為3H,3H,3H,2H,2.5H,2H.同樣L下劉晶波等提出的黏彈性邊界(VSB-Liu)的結果也被展示在圖5 中.從圖中可以看出,本文方法的結果和參考解幾乎完全吻合,而黏彈性邊界基本是在0.5 s 之后結果的精度較低.

為了研究圖3 中狄拉克脈沖的寬度對本文方法的計算精度和收斂性的影響,圖6 展示了Case-6 在T=0.2 和T=0.1 兩種載荷下,本文方法(J=1～5,L=H～5H)和VSB-Liu 方法的計算結果,圖中縱坐標為式(23)所示的A點水平加速度時程的相對誤差.T=0.2 時,本文方法中模態數n取為5;T=0.1時,n取為10.可以看出,圖6(a)(T=0.2)和圖6(b)(T=0.1)中結果的規律和圖4 (T=0.4)中Case-6結果的規律相似,都是黏彈性邊界結果遠大于本文方法的結果,本文方法J=3 時結果已收斂,再增大J值,基本不再改善精度;因近似損失的精度需要通過增大人工邊界到興趣域的距離L來彌補,都是當L=2H時,能滿足工程精度要求(相對誤差<5%),隨著L的增大,結果向0 收斂.因此,圖3 中狄拉克脈沖的寬度基本不會影響本文方法的計算精度和收斂性.

圖5 本文方法(J=3,Case-1～Case-6 對應的L 值分別為3H,3H,3H,2H,2.5H,2H)的A 點水平加速度時程結果以及同樣L 下黏彈性邊界的結果Fig.5 Time histories of horizontal acceleration at point A from the proposed method(J=3 and L=3H,3H,3H,2H,2.5H,2H for Case-1～Case-6,respectively)and from the viscous-spring boundary

圖6 T=0.2 和T=0.1 兩種狄拉克脈沖下,本文方法(J=1～5,L= H～5H)和VSB-Liu 方法的計算結果Fig.6 Results from the proposed method(J=1～5,L= H～5H)and from the VSB-Liu method under two kinds of Dirac pulse with T=0.2 and 0.1

由于本文的ABC 是基于全頻域收斂的連分式建立的,理論上ABC 應該是時域穩定的.通過后驗的方法進一步驗證本文方法的穩定性.根據線性系統理論,若ABC 系統的特征值全部分布在復平面左側,則說明ABC 穩定.圖7 展示了Case-1～Case-6 的水平方向上和豎直方向上式(19)所示的模態空間下的ABC 系統的特征值在復平面的分布.總特征值的個數為n×(J+1).從圖中可以看出ABC 特征值的實部都小于0,且其分布規律為:特征值分布在一系列半圓上,且半圓的數目是模態數n,每個半圓上特征值的數目是J+1.由于本算例每一個虛部為0 的特征值都是重根,因此圖中呈現出來的是每個半圓上特征值的數目為J.

圖7 Case-1～Case-6 的水平方向上和豎直方向上模態空間下ABC 系統的特征值在復平面的分布Fig.7 Distribution of eigenvalues in complex plane of ABC system in the horizontal and the vertical directions in the modal space for Case-1～Case-6

4 結論

本文將文獻[33]中針對標量波提出的方法擴展到矢量波,建立了一種穩定的近似時域人工邊界條件(ABC)來模擬含有瑞利阻尼的線彈性多層波導模型的平面內矢量波動.相較于文獻[33]的方法只能模擬標量波問題,本文提出的方法不僅適用于矢量波問題,對標量波問題也同樣適用.

提出的方法中影響計算精度和計算效率的參數有3 個,分別為無限域的模態數n,連分式階數J,和人工邊界遠離興趣域的距離L.數值算例表明,僅用能被載荷激起的無限域的模態數n參與計算即可.由于本文方法是一種近似方法,一般當J>3 時,再增大J值,基本不再改善精度.作者建議使用本文方法時參數J可以統一取為3.因近似處理損失的精度需要通過增大人工邊界到興趣域的距離L來彌補.確定人工邊界位置L以加速度為分析指標更為嚴格.L的取值基本與地下結構尺寸無關,它與土層的總層高H約成正比關系,關系系數與土層的材料參數有關.一般取L=3H都能將加速度時程相對誤差控制在5%以內.本文數值算例結果表明在同樣的L下,與工程中被廣泛使用的黏性邊界和黏彈性邊界相比,本文方法能夠在幾乎不增加計算量的情況下一般大約能使精度提高2～6 倍,精度的提高程度與土層的總厚度有關.

本文提出的ABC 在理論上是穩定的,數值算例也進一步驗證了其穩定性.算例表明模態空間下ABC 系統的特征值呈一定規律分布在復平面的左側,其分布規律為:特征值分布在一系列半圓上,半圓的數目是模態數n,每個半圓上特征值的數目是J+1.

本文提出的ABC 雖然人工邊界上單個結點的水平自由度和豎向自由度是空間解耦的,但是人工邊界結點間的水平自由度和結點間的豎向自由度都是空間耦聯的.下一步工作可以考慮將方法進一步優化,使人工邊界上結點間的水平自由度和結點間的豎向自由度都空間解耦,更方便工程應用.另外,可以考慮借鑒本文方法的研究思路建立矢量彈性波精確的時域人工邊界條件.