吳利華 趙 密 杜修力

(北京工業(yè)大學(xué)城市與工程安全減災(zāi)教育部重點(diǎn)實(shí)驗(yàn)室,北京 100124)

引言

用有限元法等數(shù)值方法分析大型結(jié)構(gòu)的瞬態(tài)動(dòng)力響應(yīng),一般的做法是引入人工邊界,將結(jié)構(gòu)?無(wú)限地基系統(tǒng)劃分為有限域和無(wú)限域.無(wú)限域被截去,有限域用有限元法等數(shù)值方法模擬.為了消除波在人工邊界上的虛假反射,需要在有限域的人工邊界上施加人工邊界條件(artificial boundary condition,ABC)來(lái)描述波動(dòng)在無(wú)限域中的傳播[1].ABC 是基于無(wú)限域的控制方程推導(dǎo)得到的.一般假定無(wú)限域是線性的,因?yàn)橥ǔＵJ(rèn)為無(wú)限域中只存在向無(wú)窮遠(yuǎn)傳播的波,不存在向有限域傳播的波.

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

本文的研究對(duì)象是具有確定尺寸參數(shù)的多層波導(dǎo)結(jié)構(gòu).Wang 等[4]研究表明,當(dāng)多層波導(dǎo)內(nèi)部結(jié)構(gòu)尺寸存在隨機(jī)缺陷時(shí),彈性波會(huì)產(chǎn)生局部化現(xiàn)象.該現(xiàn)象使得在原本可以傳播的頻域范圍內(nèi)彈性波被禁止傳播,進(jìn)而多層波導(dǎo)結(jié)構(gòu)中會(huì)產(chǎn)生不均勻的應(yīng)力分布.鑒于問(wèn)題的復(fù)雜性以及篇幅的限制,本文的研究范圍不包括上述多層波導(dǎo)中隨機(jī)失諧引起的彈性波局部化問(wèn)題.

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

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

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

本文的目標(biāo)是建立多層波導(dǎo)矢量波穩(wěn)定的時(shí)域ABC.借鑒文獻(xiàn)[34]的做法將多層波導(dǎo)矢量波動(dòng)方程強(qiáng)行解耦,得到x方向和y方向的兩個(gè)標(biāo)量波動(dòng)方程.將文獻(xiàn)[33]構(gòu)建達(dá)朗貝爾黏彈性多層波導(dǎo)標(biāo)量波的時(shí)域ABC 的思路用于解耦的x方向和y方向的兩個(gè)標(biāo)量波動(dòng)方程,構(gòu)建含有瑞利阻尼的線彈性多層波導(dǎo)矢量波的穩(wěn)定的時(shí)域ABC.

1 問(wèn)題描述

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

圖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,用來(lái)表示被截掉的無(wú)限域和有限域的相互作用.

一般認(rèn)為無(wú)限域是線性的,若無(wú)限域包含非線性,需要對(duì)其進(jìn)行等效線性化處理,使其滿足線性條件.無(wú)限域的波動(dòng)方程為

其中,變量上方的點(diǎn)表示對(duì)時(shí)間求導(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ù)可不同.另外,無(wú)限域的上表面物理邊界條件是應(yīng)力為0,下表面物理邊界條件是位移為0.土層的初始狀態(tài)為靜止.

2 時(shí)域人工邊界條件

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

2.1 彈性波動(dòng)方程簡(jiǎn)化

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

其中,x和y解耦簡(jiǎn)化處理后的應(yīng)力?應(yīng)變關(guān)系為

經(jīng)過(guò)解耦簡(jiǎn)化處理后的上下表面物理邊界條件分別為

其中,H為半無(wú)限域的總層高.

2.2 模態(tài)空間下的頻域動(dòng)力剛度

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

沿著土層y方向半離散簡(jiǎn)化處理后的無(wú)限域波動(dòng)方程(3),同時(shí)考慮簡(jiǎn)化后的物理邊界條件式(4)及介質(zhì)交界面的連續(xù)條件,得到無(wú)限域簡(jiǎn)化彈性波動(dòng)方程的比例邊界有限元(SBFE)方程

基于簡(jiǎn)化后的本構(gòu)關(guān)系得到無(wú)限域的人工邊界上x和y方向的等效結(jié)點(diǎn)力分別為

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

為了使問(wèn)題簡(jiǎn)單化,利用模態(tài)分解法[36],將物理空間下的SBFE 方程變換到模態(tài)空間.計(jì)算如下廣義特征值問(wèn)題

其中,上標(biāo)?1 表示對(duì)矩陣求逆.

2.3 動(dòng)力剛度的連分式展開

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

由于x方向和y方向的SBFE 方程形式相同,以下僅給出x方向的推導(dǎo)過(guò)程,y方向的推導(dǎo)同x方向.以下推導(dǎo)思路類似文獻(xiàn)[33],本文僅給出與文獻(xiàn)[33]不同的地方以及相關(guān)重要公式,詳細(xì)推導(dǎo)過(guò)程可參見(jiàn)文獻(xiàn)[33].根據(jù)式(12)和式(13),無(wú)限域x方向模態(tài)空間的動(dòng)力剛度可表示為如下的連分式

和文獻(xiàn)[33]一樣,連分式(14)也可用來(lái)近似地表示多層介質(zhì)的動(dòng)力剛度.對(duì)于多層介質(zhì),連分式中的通過(guò)求解如下黎卡提方程得到

2.4 基于連分式的時(shí)域吸收邊界條件

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

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

式(20)中,人工邊界上單個(gè)結(jié)點(diǎn)的水平自由度和豎向自由度是空間解耦的,但是人工邊界結(jié)點(diǎn)間的水平自由度和結(jié)點(diǎn)間的豎向自由度都是空間耦聯(lián)的.其中

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

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

式(21)可通過(guò)標(biāo)準(zhǔn)的時(shí)間積分算法如Newmark 常平均加速度法等求解.

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

其中,nL和nR分別表示左無(wú)限域和右無(wú)限域參與計(jì)算的模態(tài)數(shù).J表示連分式的階數(shù),通過(guò)算例分析統(tǒng)計(jì),一般可取J=3,J的取值會(huì)在第3 節(jié)詳細(xì)闡述.對(duì)于實(shí)際地震工程等問(wèn)題,根據(jù)經(jīng)驗(yàn)一般取模態(tài)數(shù)n≤3.所以額外增加的自由度數(shù)一般不會(huì)超過(guò)36,與有限域的自由度數(shù)相比,額外增加的自由度數(shù)幾乎可以忽略不計(jì).

本文建立的時(shí)域ABC 是對(duì)無(wú)限域矢量波動(dòng)的近似模擬,為了實(shí)現(xiàn)高精度,需要將有限域的人工邊界放置在距離興趣域L處,L的取值將在第3 節(jié)通過(guò)數(shù)值算例詳細(xì)分析.

本文的ABC 是基于全頻域收斂穩(wěn)定的連分式建立的,所以理論上ABC 應(yīng)該是時(shí)域穩(wěn)定的.方法的穩(wěn)定性將在第3 節(jié)通過(guò)數(shù)值算例進(jìn)一步說(shuō)明.

3 方法的性能

本節(jié)通過(guò)數(shù)值算例說(shuō)明方法的精度和穩(wěn)定性.本文方法包含3 個(gè)可變參數(shù):無(wú)限域的模態(tài)數(shù)n,連分式階數(shù)J和人工邊界遠(yuǎn)離興趣域的距離L.文獻(xiàn)[33]表明,僅用能被載荷激起的無(wú)限域的模態(tài)數(shù)n參與計(jì)算,不但能減少計(jì)算量,同時(shí)還不影響精度.下文詳細(xì)分析J和L的取值要求.

分析如圖2 所示的3 種模型.Model 1 和Model 2 都是自由場(chǎng)地,Model 1 的總層高H為20 m,Model 2 的總層高H為40 m.比較Model 1 和Model 2 的結(jié)果可以考察L,J取值與H的關(guān)系.Model 3 內(nèi)含一個(gè)尺寸為5×5 的地下方洞,其總層高H為20 m.比較Model 1 和Model 3 的結(jié)果可以考察有無(wú)地下結(jié)構(gòu)及地下結(jié)構(gòu)尺寸對(duì)L,J取值的影響.3 種模型都是在地表作用10 m 長(zhǎng)的水平線載荷,線載荷為狄拉克脈沖,其時(shí)程及頻譜如圖3 所示.為了考察狄拉克脈沖寬度是否會(huì)影響本文方法的精度和收斂性,考慮三種狄拉克脈沖:T=0.4,0.2 和0.1,其相應(yīng)的頻譜范圍分別為0～10 Hz,0～20 Hz 和0～40 Hz.圖2 中每個(gè)模型的興趣域(region of interest,ROI)為藍(lán)虛線框起來(lái)的部分,紅實(shí)線表示人工邊界,兩條紅實(shí)線框起來(lái)的部分是用于計(jì)算的有限域.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ù)對(duì)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ù)見(jiàn)表1.6 種工況的計(jì)算時(shí)長(zhǎng)都為2 s,輸出A 點(diǎn)的水平位移時(shí)程和水平加速度時(shí)程.

用有限元法分析Case-1～Case-6.圖2 所示的Model 1～Model 3 中兩條紅實(shí)線表示的人工邊界框起來(lái)的部分為有限域,人工邊界到興趣域的距離被標(biāo)記為L(zhǎng),在人工邊界處施加ABC.有限域用四邊形網(wǎng)格離散,Newmark 常平均加速度法用于時(shí)間積分計(jì)算.狄拉克脈沖中T=0.4 時(shí),有限域的網(wǎng)格尺寸為2.5 m×2.5 m,時(shí)間步長(zhǎng)為0.005 s.T=0.2 和0.1 時(shí)的網(wǎng)格尺寸和時(shí)間步長(zhǎng)分別是T=0.4 時(shí)的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ù)的頻率點(diǎn)的選擇方式有較多研究[37].本節(jié)的重點(diǎn)是驗(yàn)證本文提出方法的精度和穩(wěn)定性,由于瑞利阻尼系數(shù)的具體數(shù)值不影響本文方法的性能,所以這里采用簡(jiǎn)單的頻率點(diǎn)的選取方式,取土層前2 階頻率用于確定其瑞利阻尼系數(shù)a0和a1.Case-1～Case-6 土層的阻尼比ξ 都取為0.05,表1 列出了6 種工況的瑞利阻尼系數(shù)值.狄拉克脈沖中T=0.4 時(shí),Case-1～Case-6 用于計(jì)算的無(wú)限域模態(tài)數(shù)n分別為2,3,2,1,3 和1.

將足夠大有限域的有限元結(jié)果作為驗(yàn)證本文方法的參考解,大有限域尺寸的選取原則是人工邊界上的虛假反射波在計(jì)算時(shí)間內(nèi)不會(huì)污染興趣域.

用數(shù)值算例結(jié)果分析J和L的取值對(duì)方法精度的影響.定義如下時(shí)程相對(duì)誤差

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

圖4 展示了Case-1～Case-6 在連分式階數(shù)J分別為1,2,3,4,5 的情況下,L=H～5H時(shí)對(duì)應(yīng)的A點(diǎn)水平加速度時(shí)程的相對(duì)誤差,這里狄拉克脈沖中T=0.4.為了說(shuō)明本文方法的優(yōu)越性,圖4 也展示了同樣L下黏性邊界[5](VB)、劉晶波等[9]提出的黏彈性邊界(VSB-Liu)和杜修力等[10]提出的黏彈性邊界(VSB-Du)的結(jié)果.黏性邊界和黏彈性邊界都是經(jīng)典的ABCs,因其形式簡(jiǎn)單且時(shí)域穩(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 個(gè)公式中每層的材料參數(shù)按實(shí)際值取,r認(rèn)為是有限域的豎向?qū)ΨQ軸到人工邊界的水平距離.從圖中可以看出,6 種工況呈現(xiàn)的結(jié)果規(guī)律類似.每種工況下,黏性邊界和2 種黏彈性邊界結(jié)果的相對(duì)誤差接近,都遠(yuǎn)大于本文方法的結(jié)果.本文方法的J=1,2,3,4,5 對(duì)應(yīng)的5 條曲線都是J=1 時(shí)相對(duì)誤差較大,J=2 時(shí)相對(duì)誤差快速變小,J=3,4 和5時(shí)三者的結(jié)果基本重合.根據(jù)大量算例結(jié)果(不僅限于圖4 所列的算例結(jié)果)統(tǒng)計(jì),當(dāng)J>3 時(shí),再增加J的值,基本不再改善精度.這是因?yàn)楸疚姆椒ㄊ且环N近似方法,連分式只是近似地表示無(wú)限域的動(dòng)力剛度,即使取很高的階數(shù),連分式也不能逼近無(wú)限域動(dòng)力剛度的精確解.因近似損失的精度需要通過(guò)增大人工邊界到興趣域的距離L來(lái)彌補(bǔ).為了減少可變因素,作者建議在使用本文方法時(shí),參數(shù)J可以統(tǒng)一取為3.

圖4 6 種工況下基于本文方法、黏性邊界(VB)和2 種黏彈性邊界(VSB-Liu 和VSB-Du)的A 點(diǎn)水平加速度時(shí)程的相對(duì)誤差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

進(jìn)一步詳細(xì)分析L的選取原則.分別計(jì)算Case-1～Case-6,當(dāng)J=3,L為H～5H時(shí)A點(diǎn)水平位移時(shí)程和水平加速度時(shí)程各自的相對(duì)誤差,這里狄拉克脈沖中T=0.4.將相對(duì)誤差≤5%時(shí)各自對(duì)應(yīng)的L值,以及同樣L下黏性邊界和黏彈性邊界結(jié)果的相對(duì)誤差列于表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)對(duì)比表2 的最后3 列可以看出,黏性邊界和黏彈性邊界結(jié)果的相對(duì)誤差接近,都約為本文方法結(jié)果的5～6 倍左右,除了Case-2 和Case-5 的加速度時(shí)程結(jié)果,其黏性邊界和黏彈性邊界結(jié)果的相對(duì)誤差是本文方法結(jié)果的2～3 倍.另外,本文方法與黏性邊界和黏彈性邊界相比,僅多出少量的輔助自由度數(shù),對(duì)于Case-1～Case-6,其值分別為24,36,24,12,36,12.因此本文方法相比于黏性邊界和黏彈性邊界,能夠在幾乎不增加計(jì)算量的情況下使精度提高約2～6倍.精度的提高程度與土層的總厚度有關(guān),原因是本文方法是針對(duì)層模型建立的,其應(yīng)用在深厚土層和淺土層都具有高精度;黏性邊界和黏彈性邊界是針對(duì)半空間模型建立的,由于深厚土層相較于淺土層更接近半空間模型,因此其在深厚土層中要比在淺土層中的精度高.

(2)對(duì)比表2 第3 列表示的水平位移對(duì)應(yīng)的L值和水平加速度對(duì)應(yīng)的L值可以看出,每種工況下,水平加速度對(duì)應(yīng)的L值都是水平位移對(duì)應(yīng)結(jié)果的2 倍,除了Case-5 是2.5 倍.說(shuō)明人工邊界的位置以加速度為控制指標(biāo)來(lái)確定更為嚴(yán)格.

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

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

表2 本文方法的A 點(diǎn)水平位移時(shí)程和水平加速度時(shí)程各自的相對(duì)誤差≤5%時(shí)對(duì)應(yīng)的L 值,以及同樣L 下黏性邊界和2 種黏彈性邊界結(jié)果的相對(duì)誤差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 展示了本文方法的計(jì)算結(jié)果滿足工程精度要求(相對(duì)誤差≤5%)時(shí)的A點(diǎn)水平加速度時(shí)程.其中,J=3,Case-1～Case-6 對(duì)應(yīng)的L值分別為3H,3H,3H,2H,2.5H,2H.同樣L下劉晶波等提出的黏彈性邊界(VSB-Liu)的結(jié)果也被展示在圖5 中.從圖中可以看出,本文方法的結(jié)果和參考解幾乎完全吻合,而黏彈性邊界基本是在0.5 s 之后結(jié)果的精度較低.

為了研究圖3 中狄拉克脈沖的寬度對(duì)本文方法的計(jì)算精度和收斂性的影響,圖6 展示了Case-6 在T=0.2 和T=0.1 兩種載荷下,本文方法(J=1～5,L=H～5H)和VSB-Liu 方法的計(jì)算結(jié)果,圖中縱坐標(biāo)為式(23)所示的A點(diǎn)水平加速度時(shí)程的相對(duì)誤差.T=0.2 時(shí),本文方法中模態(tài)數(shù)n取為5;T=0.1時(shí),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 時(shí)結(jié)果已收斂,再增大J值,基本不再改善精度;因近似損失的精度需要通過(guò)增大人工邊界到興趣域的距離L來(lái)彌補(bǔ),都是當(dāng)L=2H時(shí),能滿足工程精度要求(相對(duì)誤差<5%),隨著L的增大,結(jié)果向0 收斂.因此,圖3 中狄拉克脈沖的寬度基本不會(huì)影響本文方法的計(jì)算精度和收斂性.

圖5 本文方法(J=3,Case-1～Case-6 對(duì)應(yīng)的L 值分別為3H,3H,3H,2H,2.5H,2H)的A 點(diǎn)水平加速度時(shí)程結(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 方法的計(jì)算結(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)該是時(shí)域穩(wěn)定的.通過(guò)后驗(yàn)的方法進(jìn)一步驗(yàn)證本文方法的穩(wěn)定性.根據(jù)線性系統(tǒng)理論,若ABC 系統(tǒng)的特征值全部分布在復(fù)平面左側(cè),則說(shuō)明ABC 穩(wěn)定.圖7 展示了Case-1～Case-6 的水平方向上和豎直方向上式(19)所示的模態(tài)空間下的ABC 系統(tǒng)的特征值在復(fù)平面的分布.總特征值的個(gè)數(shù)為n×(J+1).從圖中可以看出ABC 特征值的實(shí)部都小于0,且其分布規(guī)律為:特征值分布在一系列半圓上,且半圓的數(shù)目是模態(tài)數(shù)n,每個(gè)半圓上特征值的數(shù)目是J+1.由于本算例每一個(gè)虛部為0 的特征值都是重根,因此圖中呈現(xiàn)出來(lái)的是每個(gè)半圓上特征值的數(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

4 結(jié)論

本文將文獻(xiàn)[33]中針對(duì)標(biāo)量波提出的方法擴(kuò)展到矢量波,建立了一種穩(wěn)定的近似時(shí)域人工邊界條件(ABC)來(lái)模擬含有瑞利阻尼的線彈性多層波導(dǎo)模型的平面內(nèi)矢量波動(dòng).相較于文獻(xiàn)[33]的方法只能模擬標(biāo)量波問(wèn)題,本文提出的方法不僅適用于矢量波問(wèn)題,對(duì)標(biāo)量波問(wèn)題也同樣適用.

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

本文提出的ABC 在理論上是穩(wěn)定的,數(shù)值算例也進(jìn)一步驗(yàn)證了其穩(wěn)定性.算例表明模態(tài)空間下ABC 系統(tǒng)的特征值呈一定規(guī)律分布在復(fù)平面的左側(cè),其分布規(guī)律為:特征值分布在一系列半圓上,半圓的數(shù)目是模態(tài)數(shù)n,每個(gè)半圓上特征值的數(shù)目是J+1.

本文提出的ABC 雖然人工邊界上單個(gè)結(jié)點(diǎn)的水平自由度和豎向自由度是空間解耦的,但是人工邊界結(jié)點(diǎn)間的水平自由度和結(jié)點(diǎn)間的豎向自由度都是空間耦聯(lián)的.下一步工作可以考慮將方法進(jìn)一步優(yōu)化,使人工邊界上結(jié)點(diǎn)間的水平自由度和結(jié)點(diǎn)間的豎向自由度都空間解耦,更方便工程應(yīng)用.另外,可以考慮借鑒本文方法的研究思路建立矢量彈性波精確的時(shí)域人工邊界條件.