吳利華 趙 密 杜修力
(北京工業(yè)大學(xué)城市與工程安全減災(zāi)教育部重點實驗室,北京 100124)
用有限元法等數(shù)值方法分析大型結(jié)構(gòu)的瞬態(tài)動力響應(yīng),一般的做法是引入人工邊界,將結(jié)構(gòu)?無限地基系統(tǒng)劃分為有限域和無限域.無限域被截去,有限域用有限元法等數(shù)值方法模擬.為了消除波在人工邊界上的虛假反射,需要在有限域的人工邊界上施加人工邊界條件(artificial boundary condition,ABC)來描述波動在無限域中的傳播[1].ABC 是基于無限域的控制方程推導(dǎo)得到的.一般假定無限域是線性的,因為通常認(rèn)為無限域中只存在向無窮遠(yuǎn)傳播的波,不存在向有限域傳播的波.
剛性基巖上臥成層土是一種普遍存在的地基形式,又稱其為半開放的波導(dǎo)模型.與開放的半空間模型中只存在傳向無窮遠(yuǎn)的行波相比,波導(dǎo)模型由于截止頻率的存在,介質(zhì)中同時存在截止頻率以下的近場快衰波和截止頻率以上的傳向無窮遠(yuǎn)的行波[2-3].對于半空間模型,行波的振幅隨著傳播距離的增大快速減弱,所以即使不考慮介質(zhì)材料阻尼的吸收效應(yīng),其無限域也可以被假定為小應(yīng)變,即可以被假定為線性的.而對于波導(dǎo)模型,由于能量被困在介質(zhì)中,行波的振幅衰減得較為緩慢,如果僅考慮幾何散射不考慮介質(zhì)吸收效應(yīng),將無限域假定為小應(yīng)變是有問題的,所以更合理的做法應(yīng)該是將波導(dǎo)模型考慮為有材料阻尼的線彈性介質(zhì).
本文的研究對象是具有確定尺寸參數(shù)的多層波導(dǎo)結(jié)構(gòu).Wang 等[4]研究表明,當(dāng)多層波導(dǎo)內(nèi)部結(jié)構(gòu)尺寸存在隨機缺陷時,彈性波會產(chǎn)生局部化現(xiàn)象.該現(xiàn)象使得在原本可以傳播的頻域范圍內(nèi)彈性波被禁止傳播,進而多層波導(dǎo)結(jié)構(gòu)中會產(chǎn)生不均勻的應(yīng)力分布.鑒于問題的復(fù)雜性以及篇幅的限制,本文的研究范圍不包括上述多層波導(dǎo)中隨機失諧引起的彈性波局部化問題.
如有限域含有非線性等許多情況下,需要ABC是時域方法.相較于半空間模型的時域ABC 已有大量的研究工作[5-13],波導(dǎo)模型的時域ABC 研究工作較少.由于波導(dǎo)模型的ABC 需要能夠同時模擬快衰波和行波的能量傳遞,直接在時域下建立其ABC 是有挑戰(zhàn)的.根據(jù)是否精確地處理波導(dǎo)模型無限域定解問題,以下將已有的波導(dǎo)模型的ABC 分為精確的ABC 和近似的ABC 進行闡述.
精確的ABC 一般是首先基于解析法或半解析半離散法獲得精確的頻域動力剛度,再通過時域化方法建立時域ABC.解析法一般對單層介質(zhì)標(biāo)量波[14-16]適用,多層介質(zhì)的解析解難以獲得.半解析半離散法只對無限域的人工邊界進行離散,對于單層介質(zhì)和多層介質(zhì)都適用,如邊界元法(BEM)[17]、薄層法(TLM)[18-21]和比例邊界有限元法(SBFEM)[3,22-26].BEM 需要基本解,而多層介質(zhì)格林函數(shù)的基本解是復(fù)雜的,基本解的存在與否限制了BEM 在工程實際中的應(yīng)用.相比于BEM,TLM 和SBFEM 不需要基本解,能夠自動滿足無窮遠(yuǎn)輻射邊界條件.將頻域動力剛度變換到時域,一種方法是時間卷積計算[3,25],該方法不僅耗時,而且不適用于非線性分析.學(xué)者們通過時間局部化方法建立了精確穩(wěn)定的時間局部的ABC.Zhao 等[15]基于解析法獲得2D 單層波導(dǎo)標(biāo)量波的動力剛度,再用有理函數(shù)逼近方法建立精確的時域ABC.Prempramote 等[22]將解析法獲得的動力剛度用雙漸近連分式展開,建立了2D 單層標(biāo)量波的時域ABC.Li 等[16]基于解析法獲得2D 和3D 均勻?qū)訕?biāo)量波的動力剛度,用動力剛度連分式展開的時間局部化方法建立了精確的時域ABC.Liu 等[26]基于半解析半離散法,利用算子分離法獲得2D 單層標(biāo)量波的時域ABC.Wang 等[23]基于SBFEM 和動力剛度的雙漸近連分式展開[22]建立了2D 庫水的時域ABC,并將方法嵌入到通用有限元軟件ABAQUS 中用來計算壩?庫水的動力相互作用.但是上述精確的時域ABC 都是針對單層波導(dǎo)標(biāo)量波獲得.根據(jù)作者的調(diào)查,目前幾乎沒有關(guān)于多層波導(dǎo)穩(wěn)定的精確時域ABC 的相關(guān)報導(dǎo).因為對于多層波導(dǎo),其精確的時域ABC 易發(fā)生失穩(wěn),目前還沒有較好地抑制其失穩(wěn)的方法.
多層波導(dǎo)精確的時域ABC 難以獲得,學(xué)者們通過對無限域定解問題作適量的簡化,建立了多層波導(dǎo)近似的時域ABC.近似的低階時空局部的ABC 由于其形式簡單,在工程中被廣泛使用,如黏性邊界[5,27]、多次透射邊界[7,11]以及黏彈性邊界[9-10].這類低階時空局部ABC 是基于半空間模型建立的,用于波導(dǎo)模型時精度較低.Hagstrom 等[28]建立了多層波導(dǎo)標(biāo)量波的高階時空局部ABC.完美匹配層法[29-30]是近年的研究熱點,它是通過在人工邊界外附加一層具有耗能作用的緩沖區(qū)來吸收透射波,但吸收層具有一定厚度,層厚度和耗能參數(shù)選擇不當(dāng)時,會導(dǎo)致精度低或發(fā)生失穩(wěn).Zhao 等[32]基于半解析半離散法和一種全頻域收斂的動力剛度連分式展開構(gòu)建了多層波導(dǎo)標(biāo)量波的時域ABC,并對其穩(wěn)定性進行了證明.吳利華等[33]將多層波導(dǎo)考慮為達(dá)朗貝爾黏彈性介質(zhì),利用和文獻[32]類似的思路建立了考慮阻尼的多層波導(dǎo)標(biāo)量波的時域ABC.高毅超等[34]將彈性動力學(xué)矢量方程強行解耦為兩個標(biāo)量波方程,基于SBFEM和雙漸近連分式[22]構(gòu)建了一種模擬單層波導(dǎo)矢量波的近似ABC.Liu 等[35]將文獻[34]的方法嵌入到開源有限元軟件OpenSees 中,分析了地下車站的地震動響應(yīng).Liu 等[31]基于半解析半離散法和算子分離法建立了多層波導(dǎo)矢量波的近似ABC,但該方法有時會出現(xiàn)失穩(wěn)現(xiàn)象.就作者所知,目前幾乎沒有專門針對多層波導(dǎo)矢量波的穩(wěn)定時域ABC 的相關(guān)研究.
本文的目標(biāo)是建立多層波導(dǎo)矢量波穩(wěn)定的時域ABC.借鑒文獻[34]的做法將多層波導(dǎo)矢量波動方程強行解耦,得到x方向和y方向的兩個標(biāo)量波動方程.將文獻[33]構(gòu)建達(dá)朗貝爾黏彈性多層波導(dǎo)標(biāo)量波的時域ABC 的思路用于解耦的x方向和y方向的兩個標(biāo)量波動方程,構(gòu)建含有瑞利阻尼的線彈性多層波導(dǎo)矢量波的穩(wěn)定的時域ABC.
用有限元法分析如圖1 所示的剛性基巖上臥成層地基上結(jié)構(gòu)在動力載荷作用下平面內(nèi)的時程響應(yīng).由于有限元法只能分析有限區(qū)域,需要引入人工邊界將左右兩個半無限域截去.為了避免波在人工邊界處引起虛假反射,需要在有限域的人工邊界處施加一個人工邊界條件(ABC)來模擬波在無限域中的傳播.
圖1 剛性基巖上臥成層土Fig.1 A multilayered medium overlying rigid bedrock
有限域可以是非線性的,其有限元方程為
其中,下標(biāo)b 和i 分別表示人工邊界部分和有限域除了人工邊界的其他部分,M是質(zhì)量矩陣,C是阻尼矩陣,K是剛度矩陣,u是平面內(nèi)位移向量,是非線性內(nèi)力,b是體力,fi是施加在有限域的載荷,fb是施加在有限域人工邊界上的力,即本文要建立的ABC,用來表示被截掉的無限域和有限域的相互作用.
一般認(rèn)為無限域是線性的,若無限域包含非線性,需要對其進行等效線性化處理,使其滿足線性條件.無限域的波動方程為
其中,變量上方的點表示對時間求導(dǎo),(x,y)表示笛卡爾坐標(biāo),上標(biāo)T 表示轉(zhuǎn)置;U={ux uy}T,ux,uy分別為x方向和y方向的位移,ux和uy在不同土層的分界面以及人工邊界處是連續(xù)的;ρ 是質(zhì)量密度,λ 是拉梅常數(shù),G是剪切模量,每一層的材料參數(shù)是常量,不同層的材料參數(shù)可不同.另外,無限域的上表面物理邊界條件是應(yīng)力為0,下表面物理邊界條件是位移為0.土層的初始狀態(tài)為靜止.
將無限域作為研究對象,建立時域吸收邊界條件.文獻[33]建立了模擬剛性基巖上臥多層介質(zhì)中標(biāo)量波傳播的ABC,本文將該方法擴展到矢量波.方法的思路是:(1)將矢量波動方程簡化為x方向和y方向解耦的兩個標(biāo)量波動方程;(2)基于比例邊界有限元法,得到兩個標(biāo)量波動方程模態(tài)空間下半離散的頻域動力剛度方程;(3)通過文獻[33]中提出的連分式來表示頻域動力剛度;(4)引入輔助變量,將連分式時域化,得到時域吸收邊界條件;(5)所建立的ABC 可以和有限域的有限元方程無縫耦合,通過數(shù)值積分算法求解.
高毅超等[34]通過將矢量波動方程強行解耦,基于雙漸近連分式建立了單層波導(dǎo)模型的ABC.本文借鑒其做法,忽略無限域波動方程(2)中x和y的耦合項,簡化后得到關(guān)于x方向和y方向的兩個標(biāo)量方程分別為
其中,x和y解耦簡化處理后的應(yīng)力?應(yīng)變關(guān)系為
經(jīng)過解耦簡化處理后的上下表面物理邊界條件分別為
其中,H為半無限域的總層高.
根據(jù)比例邊界有限元法(SBFEM)[23]建立模態(tài)空間下的頻域動力剛度方程.
沿著土層y方向半離散簡化處理后的無限域波動方程(3),同時考慮簡化后的物理邊界條件式(4)及介質(zhì)交界面的連續(xù)條件,得到無限域簡化彈性波動方程的比例邊界有限元(SBFE)方程
基于簡化后的本構(gòu)關(guān)系得到無限域的人工邊界上x和y方向的等效結(jié)點力分別為
引言的第二段闡述了波導(dǎo)模型應(yīng)該考慮阻尼的吸收效應(yīng).由于瑞利阻尼是目前應(yīng)用最廣泛的阻尼矩陣的數(shù)值模型,本文假定無限域土層含有瑞利阻尼.瑞利阻尼定義阻尼矩陣是質(zhì)量陣和剛度陣的線性組合[36],即CRayleigh=a0M+a1K,其中,ξ 為陣型阻尼比,ωi,ωj為選擇用于確定阻尼系數(shù)a0和a1的頻率點.瑞利阻尼的引入體現(xiàn)在式(5a)和式(5b)的第四項,需要說明的是,這里無限域瑞利阻尼的質(zhì)量陣和剛度陣都是通過半離散x和y解耦的標(biāo)量控制方程得到的,注意其與二維有限元法中的瑞利阻尼相區(qū)別.
為了使問題簡單化,利用模態(tài)分解法[36],將物理空間下的SBFE 方程變換到模態(tài)空間.計算如下廣義特征值問題
其中,上標(biāo)?1 表示對矩陣求逆.
其中,z=iω/ωn,β=a0/ωn+a1ωn,i 是虛數(shù)單位,ω 表示圓頻率.x和y方向分別有c=cp和c=cs.
由于x方向和y方向的SBFE 方程形式相同,以下僅給出x方向的推導(dǎo)過程,y方向的推導(dǎo)同x方向.以下推導(dǎo)思路類似文獻[33],本文僅給出與文獻[33]不同的地方以及相關(guān)重要公式,詳細(xì)推導(dǎo)過程可參見文獻[33].根據(jù)式(12)和式(13),無限域x方向模態(tài)空間的動力剛度可表示為如下的連分式
和文獻[33]一樣,連分式(14)也可用來近似地表示多層介質(zhì)的動力剛度.對于多層介質(zhì),連分式中的通過求解如下黎卡提方程得到
同文獻[33]一樣,通過引入輔助變量的方式將連分式時域化.將式(14)代入式(10a),同時引入輔助變量,得到
如今,微課的設(shè)計絕大多數(shù)都是對教材中知識點的教授,還滯留在對學(xué)生采取知識灌溉的方式上,然而微課的作用絕不止在于此。微課的時間很短,但是使用微課的最終目標(biāo)是希望做到事半功倍,以小見大。所以,教師必須從簡潔的授課型微課入手,漸漸向啟發(fā)型、探究型的教學(xué)模式發(fā)展。想要做到這一點,教師在使用微課之前,就要做好充實的資料收集與課前準(zhǔn)備,以便能夠科學(xué)有效地權(quán)衡探究型活動和師生互動使用的微課時間。
式(20)中,人工邊界上單個結(jié)點的水平自由度和豎向自由度是空間解耦的,但是人工邊界結(jié)點間的水平自由度和結(jié)點間的豎向自由度都是空間耦聯(lián)的.其中
本文建立的ABC 可與有限元法無縫耦合.將式(20)代入有限元方程式(1),得
式(21)可通過標(biāo)準(zhǔn)的時間積分算法如Newmark 常平均加速度法等求解.
由于在有限域的人工邊界上施加ABC 增加的額外自由度為
其中,nL和nR分別表示左無限域和右無限域參與計算的模態(tài)數(shù).J表示連分式的階數(shù),通過算例分析統(tǒng)計,一般可取J=3,J的取值會在第3 節(jié)詳細(xì)闡述.對于實際地震工程等問題,根據(jù)經(jīng)驗一般取模態(tài)數(shù)n≤3.所以額外增加的自由度數(shù)一般不會超過36,與有限域的自由度數(shù)相比,額外增加的自由度數(shù)幾乎可以忽略不計.
本文建立的時域ABC 是對無限域矢量波動的近似模擬,為了實現(xiàn)高精度,需要將有限域的人工邊界放置在距離興趣域L處,L的取值將在第3 節(jié)通過數(shù)值算例詳細(xì)分析.
本文的ABC 是基于全頻域收斂穩(wěn)定的連分式建立的,所以理論上ABC 應(yīng)該是時域穩(wěn)定的.方法的穩(wěn)定性將在第3 節(jié)通過數(shù)值算例進一步說明.
本節(jié)通過數(shù)值算例說明方法的精度和穩(wěn)定性.本文方法包含3 個可變參數(shù):無限域的模態(tài)數(shù)n,連分式階數(shù)J和人工邊界遠(yuǎn)離興趣域的距離L.文獻[33]表明,僅用能被載荷激起的無限域的模態(tài)數(shù)n參與計算,不但能減少計算量,同時還不影響精度.下文詳細(xì)分析J和L的取值要求.
分析如圖2 所示的3 種模型.Model 1 和Model 2 都是自由場地,Model 1 的總層高H為20 m,Model 2 的總層高H為40 m.比較Model 1 和Model 2 的結(jié)果可以考察L,J取值與H的關(guān)系.Model 3 內(nèi)含一個尺寸為5×5 的地下方洞,其總層高H為20 m.比較Model 1 和Model 3 的結(jié)果可以考察有無地下結(jié)構(gòu)及地下結(jié)構(gòu)尺寸對L,J取值的影響.3 種模型都是在地表作用10 m 長的水平線載荷,線載荷為狄拉克脈沖,其時程及頻譜如圖3 所示.為了考察狄拉克脈沖寬度是否會影響本文方法的精度和收斂性,考慮三種狄拉克脈沖:T=0.4,0.2 和0.1,其相應(yīng)的頻譜范圍分別為0~10 Hz,0~20 Hz 和0~40 Hz.圖2 中每個模型的興趣域(region of interest,ROI)為藍(lán)虛線框起來的部分,紅實線表示人工邊界,兩條紅實線框起來的部分是用于計算的有限域.L表示人工邊界到興趣域的距離.
圖2 數(shù)值算例模型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)
為了分析土層材料參數(shù)對L,J取值的影響,將圖2 中3 種模型分別都考慮為均勻介質(zhì)和分別都考慮為4 層介質(zhì).分別將單層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 土層都被考慮為有瑞利阻尼的線彈性介質(zhì),其土層材料參數(shù)見表1.6 種工況的計算時長都為2 s,輸出A 點的水平位移時程和水平加速度時程.
用有限元法分析Case-1~Case-6.圖2 所示的Model 1~Model 3 中兩條紅實線表示的人工邊界框起來的部分為有限域,人工邊界到興趣域的距離被標(biāo)記為L,在人工邊界處施加ABC.有限域用四邊形網(wǎng)格離散,Newmark 常平均加速度法用于時間積分計算.狄拉克脈沖中T=0.4 時,有限域的網(wǎng)格尺寸為2.5 m×2.5 m,時間步長為0.005 s.T=0.2 和0.1 時的網(wǎng)格尺寸和時間步長分別是T=0.4 時的0.5 倍和0.25 倍.
表1 Case-1~Case-6 的土層材料參數(shù)及其瑞利阻尼系數(shù)Table 1 Material parameters and Rayleigh damping coefficients of Case-1~Case-6
用于確定瑞利阻尼系數(shù)的頻率點的選擇方式有較多研究[37].本節(jié)的重點是驗證本文提出方法的精度和穩(wěn)定性,由于瑞利阻尼系數(shù)的具體數(shù)值不影響本文方法的性能,所以這里采用簡單的頻率點的選取方式,取土層前2 階頻率用于確定其瑞利阻尼系數(shù)a0和a1.Case-1~Case-6 土層的阻尼比ξ 都取為0.05,表1 列出了6 種工況的瑞利阻尼系數(shù)值.狄拉克脈沖中T=0.4 時,Case-1~Case-6 用于計算的無限域模態(tài)數(shù)n分別為2,3,2,1,3 和1.
將足夠大有限域的有限元結(jié)果作為驗證本文方法的參考解,大有限域尺寸的選取原則是人工邊界上的虛假反射波在計算時間內(nèi)不會污染興趣域.
用數(shù)值算例結(jié)果分析J和L的取值對方法精度的影響.定義如下時程相對誤差
其中,Xj為基于有限元法和ABC 得到的小有限域的時程結(jié)果,Yj為僅基于有限元法得到的大有限域的時程結(jié)果(即參考解),N為總時步.
圖4 展示了Case-1~Case-6 在連分式階數(shù)J分別為1,2,3,4,5 的情況下,L=H~5H時對應(yīng)的A點水平加速度時程的相對誤差,這里狄拉克脈沖中T=0.4.為了說明本文方法的優(yōu)越性,圖4 也展示了同樣L下黏性邊界[5](VB)、劉晶波等[9]提出的黏彈性邊界(VSB-Liu)和杜修力等[10]提出的黏彈性邊界(VSB-Du)的結(jié)果.黏性邊界和黏彈性邊界都是經(jīng)典的ABCs,因其形式簡單且時域穩(wěn)定,在工程中被廣泛使用.黏性邊界[5]是基于平面波假定建立的阻尼器邊界條件,其每層的阻尼元件法向參數(shù)為CN=ρcp,切向參數(shù)為CT=ρcs.黏彈性邊界是基于柱面波假定建立的并聯(lián)彈簧?阻尼器邊界條件.劉晶波等[9]提出的黏彈性邊界每層的彈簧元件法向參數(shù)為KN=2G/r,切向參數(shù)為KT=3G/(2r),每層的阻尼元件法向參數(shù)為CN=ρcp,切向參數(shù)為CT=ρcs;杜修力等[10]提出的黏彈性邊界每層的彈簧元件法向參數(shù)為KN=(λ+2G)/(3.6r),切向參數(shù)為KT=G/(3.6r),每層的阻尼元件法向參數(shù)為CN=1.1ρcp,切向參數(shù)為CT=1.1ρcs.3 個公式中每層的材料參數(shù)按實際值取,r認(rèn)為是有限域的豎向?qū)ΨQ軸到人工邊界的水平距離.從圖中可以看出,6 種工況呈現(xiàn)的結(jié)果規(guī)律類似.每種工況下,黏性邊界和2 種黏彈性邊界結(jié)果的相對誤差接近,都遠(yuǎn)大于本文方法的結(jié)果.本文方法的J=1,2,3,4,5 對應(yīng)的5 條曲線都是J=1 時相對誤差較大,J=2 時相對誤差快速變小,J=3,4 和5時三者的結(jié)果基本重合.根據(jù)大量算例結(jié)果(不僅限于圖4 所列的算例結(jié)果)統(tǒng)計,當(dāng)J>3 時,再增加J的值,基本不再改善精度.這是因為本文方法是一種近似方法,連分式只是近似地表示無限域的動力剛度,即使取很高的階數(shù),連分式也不能逼近無限域動力剛度的精確解.因近似損失的精度需要通過增大人工邊界到興趣域的距離L來彌補.為了減少可變因素,作者建議在使用本文方法時,參數(shù)J可以統(tǒng)一取為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
進一步詳細(xì)分析L的選取原則.分別計算Case-1~Case-6,當(dāng)J=3,L為H~5H時A點水平位移時程和水平加速度時程各自的相對誤差,這里狄拉克脈沖中T=0.4.將相對誤差≤5%時各自對應(yīng)的L值,以及同樣L下黏性邊界和黏彈性邊界結(jié)果的相對誤差列于表2.根據(jù)表2 可以得出以下結(jié)論(表中第2 列u表示位移,a表示加速度;第3 列表頭L表示人工邊界到興趣域的距離;第4 列表頭ABC 表示本文方法的結(jié)果;第5 列表頭VB 表示黏性邊界[5]的結(jié)果;第6 列表頭VSB-Liu 表示劉晶波等[9]提出的黏彈性邊界的結(jié)果;第7 列表頭VSB-Du 表示杜修力等[10]提出的黏彈性邊界的結(jié)果).
(1)對比表2 的最后3 列可以看出,黏性邊界和黏彈性邊界結(jié)果的相對誤差接近,都約為本文方法結(jié)果的5~6 倍左右,除了Case-2 和Case-5 的加速度時程結(jié)果,其黏性邊界和黏彈性邊界結(jié)果的相對誤差是本文方法結(jié)果的2~3 倍.另外,本文方法與黏性邊界和黏彈性邊界相比,僅多出少量的輔助自由度數(shù),對于Case-1~Case-6,其值分別為24,36,24,12,36,12.因此本文方法相比于黏性邊界和黏彈性邊界,能夠在幾乎不增加計算量的情況下使精度提高約2~6倍.精度的提高程度與土層的總厚度有關(guān),原因是本文方法是針對層模型建立的,其應(yīng)用在深厚土層和淺土層都具有高精度;黏性邊界和黏彈性邊界是針對半空間模型建立的,由于深厚土層相較于淺土層更接近半空間模型,因此其在深厚土層中要比在淺土層中的精度高.
(2)對比表2 第3 列表示的水平位移對應(yīng)的L值和水平加速度對應(yīng)的L值可以看出,每種工況下,水平加速度對應(yīng)的L值都是水平位移對應(yīng)結(jié)果的2 倍,除了Case-5 是2.5 倍.說明人工邊界的位置以加速度為控制指標(biāo)來確定更為嚴(yán)格.
(3)比較單層的Case-1 和Case-3 水平加速度對應(yīng)的L值.Case-1 和Case-3 總層高相同,Case-1 為自由場地,Case-3 含有5 m×5 m 的地下結(jié)構(gòu).兩者都是L=3H時,水平加速度時程的相對誤差能控制在5%以內(nèi).說明L的取值基本不受是否含有地下結(jié)構(gòu)或結(jié)構(gòu)尺寸影響.比較4 層的Case-4 和Case-6 能得出同樣的結(jié)論.
(4)比較單層Case-1 和Case-2 水平加速度對應(yīng)的L值.Case-1 和Case-2 都是自由場地,Case-2 的總層高H是Case-1 的2 倍.兩者都是L=3H時,水平加速度時程的相對誤差能控制在5%以內(nèi).同樣,比較4 層的Case-4 和Case-5.為了將水平加速度時程的相對誤差控制在5%以內(nèi),Case-4 要求L=2H,Case-5要求L=2.5H.說明L的取值大約與土層總層高H成正比關(guān)系,關(guān)系系數(shù)與土層的材料參數(shù)有關(guān).根據(jù)大量算例統(tǒng)計,一般L=3H都能將加速度時程相對誤差控制在5%以內(nèi).
表2 本文方法的A 點水平位移時程和水平加速度時程各自的相對誤差≤5%時對應(yīng)的L 值,以及同樣L 下黏性邊界和2 種黏彈性邊界結(jié)果的相對誤差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 展示了本文方法的計算結(jié)果滿足工程精度要求(相對誤差≤5%)時的A點水平加速度時程.其中,J=3,Case-1~Case-6 對應(yīng)的L值分別為3H,3H,3H,2H,2.5H,2H.同樣L下劉晶波等提出的黏彈性邊界(VSB-Liu)的結(jié)果也被展示在圖5 中.從圖中可以看出,本文方法的結(jié)果和參考解幾乎完全吻合,而黏彈性邊界基本是在0.5 s 之后結(jié)果的精度較低.
為了研究圖3 中狄拉克脈沖的寬度對本文方法的計算精度和收斂性的影響,圖6 展示了Case-6 在T=0.2 和T=0.1 兩種載荷下,本文方法(J=1~5,L=H~5H)和VSB-Liu 方法的計算結(jié)果,圖中縱坐標(biāo)為式(23)所示的A點水平加速度時程的相對誤差.T=0.2 時,本文方法中模態(tài)數(shù)n取為5;T=0.1時,n取為10.可以看出,圖6(a)(T=0.2)和圖6(b)(T=0.1)中結(jié)果的規(guī)律和圖4 (T=0.4)中Case-6結(jié)果的規(guī)律相似,都是黏彈性邊界結(jié)果遠(yuǎn)大于本文方法的結(jié)果,本文方法J=3 時結(jié)果已收斂,再增大J值,基本不再改善精度;因近似損失的精度需要通過增大人工邊界到興趣域的距離L來彌補,都是當(dāng)L=2H時,能滿足工程精度要求(相對誤差<5%),隨著L的增大,結(jié)果向0 收斂.因此,圖3 中狄拉克脈沖的寬度基本不會影響本文方法的計算精度和收斂性.
圖5 本文方法(J=3,Case-1~Case-6 對應(yīng)的L 值分別為3H,3H,3H,2H,2.5H,2H)的A 點水平加速度時程結(jié)果以及同樣L 下黏彈性邊界的結(jié)果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 方法的計算結(jié)果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 應(yīng)該是時域穩(wěn)定的.通過后驗的方法進一步驗證本文方法的穩(wěn)定性.根據(jù)線性系統(tǒng)理論,若ABC 系統(tǒng)的特征值全部分布在復(fù)平面左側(cè),則說明ABC 穩(wěn)定.圖7 展示了Case-1~Case-6 的水平方向上和豎直方向上式(19)所示的模態(tài)空間下的ABC 系統(tǒng)的特征值在復(fù)平面的分布.總特征值的個數(shù)為n×(J+1).從圖中可以看出ABC 特征值的實部都小于0,且其分布規(guī)律為:特征值分布在一系列半圓上,且半圓的數(shù)目是模態(tài)數(shù)n,每個半圓上特征值的數(shù)目是J+1.由于本算例每一個虛部為0 的特征值都是重根,因此圖中呈現(xiàn)出來的是每個半圓上特征值的數(shù)目為J.
圖7 Case-1~Case-6 的水平方向上和豎直方向上模態(tài)空間下ABC 系統(tǒng)的特征值在復(fù)平面的分布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
本文將文獻[33]中針對標(biāo)量波提出的方法擴展到矢量波,建立了一種穩(wěn)定的近似時域人工邊界條件(ABC)來模擬含有瑞利阻尼的線彈性多層波導(dǎo)模型的平面內(nèi)矢量波動.相較于文獻[33]的方法只能模擬標(biāo)量波問題,本文提出的方法不僅適用于矢量波問題,對標(biāo)量波問題也同樣適用.
提出的方法中影響計算精度和計算效率的參數(shù)有3 個,分別為無限域的模態(tài)數(shù)n,連分式階數(shù)J,和人工邊界遠(yuǎn)離興趣域的距離L.數(shù)值算例表明,僅用能被載荷激起的無限域的模態(tài)數(shù)n參與計算即可.由于本文方法是一種近似方法,一般當(dāng)J>3 時,再增大J值,基本不再改善精度.作者建議使用本文方法時參數(shù)J可以統(tǒng)一取為3.因近似處理損失的精度需要通過增大人工邊界到興趣域的距離L來彌補.確定人工邊界位置L以加速度為分析指標(biāo)更為嚴(yán)格.L的取值基本與地下結(jié)構(gòu)尺寸無關(guān),它與土層的總層高H約成正比關(guān)系,關(guān)系系數(shù)與土層的材料參數(shù)有關(guān).一般取L=3H都能將加速度時程相對誤差控制在5%以內(nèi).本文數(shù)值算例結(jié)果表明在同樣的L下,與工程中被廣泛使用的黏性邊界和黏彈性邊界相比,本文方法能夠在幾乎不增加計算量的情況下一般大約能使精度提高2~6 倍,精度的提高程度與土層的總厚度有關(guān).
本文提出的ABC 在理論上是穩(wěn)定的,數(shù)值算例也進一步驗證了其穩(wěn)定性.算例表明模態(tài)空間下ABC 系統(tǒng)的特征值呈一定規(guī)律分布在復(fù)平面的左側(cè),其分布規(guī)律為:特征值分布在一系列半圓上,半圓的數(shù)目是模態(tài)數(shù)n,每個半圓上特征值的數(shù)目是J+1.
本文提出的ABC 雖然人工邊界上單個結(jié)點的水平自由度和豎向自由度是空間解耦的,但是人工邊界結(jié)點間的水平自由度和結(jié)點間的豎向自由度都是空間耦聯(lián)的.下一步工作可以考慮將方法進一步優(yōu)化,使人工邊界上結(jié)點間的水平自由度和結(jié)點間的豎向自由度都空間解耦,更方便工程應(yīng)用.另外,可以考慮借鑒本文方法的研究思路建立矢量彈性波精確的時域人工邊界條件.