粘彈性人工邊界單元及地震動輸入方法比較研究

王 展,景立平,3,陸新宇,齊文浩

(1.中國地震局工程力學(xué)研究所 中國地震局地震工程與工程振動重點(diǎn)實驗室,黑龍江 哈爾濱 150080;2.地震災(zāi)害防治應(yīng)急管理部重點(diǎn)實驗室,黑龍江 哈爾濱 150080;3.防災(zāi)科技學(xué)院,河北 三河 065201)

0 引言

近場波動一般是復(fù)雜介質(zhì)中三維非線性問題,通常采用數(shù)值計算方法求解。發(fā)展近場波動數(shù)值模擬技術(shù)的目標(biāo)是在保證數(shù)值精度和數(shù)值穩(wěn)定性的條件下不斷提高計算效率。作為近場波動問題的重要組成部分,人工邊界理論及技術(shù)是實現(xiàn)這一目標(biāo)的必經(jīng)之路。人工邊界需要模擬無限域?qū)τ嬎銋^(qū)域的作用,使兩者對計算區(qū)域的作用一致,保證散射波透過人工邊界或者被人工邊界吸收而不會反射回計算區(qū)域。常用的人工邊界主要有：粘性邊界、透射邊界[4]、無限元邊界[5]、粘彈性邊界和自由度綁定邊界[7,24]等。DEEKS等[6]1994年推導(dǎo)了柱坐標(biāo)系下剪切波和膨脹波波動方程,建立了平面應(yīng)變邊界條件。隨后,劉晶波等[3]給出了粘彈性邊界的一種外源荷載輸入方法[13],建立了三維粘彈性人工邊界。趙密[22]通過遠(yuǎn)場平面波的格林函數(shù)解,考慮無衰波和衰減波混合并加以推廣,得出了粘彈性邊界的廣義表達(dá)形式。劉晶波等[1,23]提出了一致粘彈性邊界的概念,并在此基礎(chǔ)上提出了一致粘彈性邊界單元,這種邊界單元和粘彈性邊界有相同的模擬精度。劉晶波等[2]2018年提出了一種基于子結(jié)構(gòu)的外源地震動輸入方法,彌補(bǔ)了粘彈性邊界單元無法有效輸入外源荷載的問題。

數(shù)值模擬技術(shù)發(fā)展的同時,還應(yīng)考慮到使用者的應(yīng)用方便。粘彈性因為其明顯的優(yōu)點(diǎn)被很多人接受,但是粘彈性邊界在實際操作中需要在每個邊界節(jié)點(diǎn)上施加兩個(三維為三個)相互垂直并聯(lián)的彈簧阻尼元件,并且彈簧阻尼元件的參數(shù)與土層材料的參數(shù)和節(jié)點(diǎn)對應(yīng)的面積有關(guān),這意味著對于非均網(wǎng)格計算模型,需要逐個單元進(jìn)行復(fù)雜的大量計算,不利于推廣使用。同時,地震動輸入方法是能否得到合理數(shù)值解答的關(guān)鍵?，F(xiàn)在主流的地震動輸入方式是波動法[8],雖然模擬精度較高,但是需要在每一時步計算每個邊界節(jié)點(diǎn)的集中力,并判斷集中力的作用方向,操作復(fù)雜且繁瑣。

本文引用了一致粘彈性邊界的概念,對與粘彈性邊界有相同計算精度的粘彈性邊界單元進(jìn)行了推導(dǎo),這種邊界只需要將單元向外延伸一定寬度,并且賦予其和相鄰內(nèi)部單元材料相關(guān)的屬性即可。同時采用劉晶波提出的一種基于局部人工邊界的子結(jié)構(gòu)的地震動輸入方法,避免繁瑣的邊界節(jié)點(diǎn)力的計算和方向判斷,且適用于各種應(yīng)力型局部人工邊界。對于外延寬度的取值,李述濤等[11]對于粘彈性邊界單元計算穩(wěn)定性做了討論,給出了側(cè)邊子系統(tǒng)和角點(diǎn)子系統(tǒng)的穩(wěn)定條件,可通過此穩(wěn)定條件確定小變形下的邊界單元寬度。基于此,本文考慮大變形時粘彈性邊界單元的適用性,以底邊界輸入位移幅值umax和邊界單元寬度b的比值ε=umax/b為參數(shù)考慮此影響。數(shù)值算例表明：當(dāng)ε≤0.1時,計算結(jié)果誤差較小,并以脈沖輸入的散射場作為算例進(jìn)行了驗證。

1 一致粘彈性人工邊界

1.1 粘彈性邊界理論

粘彈性邊界是在人工邊界上施加并聯(lián)彈簧元件和阻尼元件,用來模擬無限域?qū)τ邢抻嬎銋^(qū)域的作用。其中彈簧元件和阻尼元件的參數(shù)是根據(jù)空間柱面波或球面波理論解確定。

粘彈性人工邊界的彈簧系數(shù)KB和阻尼系數(shù)CB分別為：

(1)

(2)

式中：G表示介質(zhì)的剪切模量,R表示加載點(diǎn)到人工邊界的距離,ρ表示介質(zhì)的質(zhì)量密度,Cs和Cp分別表示介質(zhì)的剪切波速和壓縮波速,αT和αN分別表示切向和法向彈簧的修正系數(shù),其取值參考劉晶波推薦值,見表1。

表1 彈簧修正系數(shù)取值Table 1 Spring correction factor

1.2 二維一致粘彈性邊界

用有限單元法對模型進(jìn)行空間離散化時,粘彈性人工邊界也隨之離散,如果將組成粘彈性邊界的彈簧阻尼元件看作一個單元,其形函數(shù)若與內(nèi)部單元相同,認(rèn)為這樣離散后的邊界為一致粘彈性邊界。

為了推導(dǎo)一致粘彈性邊界的剛度矩陣和阻尼矩陣,將內(nèi)部單元與粘彈性邊界相鄰的一條邊視作桿件單元,i和j表示單元的兩個節(jié)點(diǎn),2L表示桿長。此時,認(rèn)為最外側(cè)桿單元的形函數(shù)和由彈簧阻尼元件組成邊界的形函數(shù)相同。

(3)

假設(shè)法向和切向有相同的形函數(shù),桿件單元的形函數(shù)矩陣為：

(4)

邊界上分布的彈簧的彈性剛度矩陣為：

(5)

根據(jù)虛功方程δW外=δW變,注意到不同于一般的單元,連續(xù)分布彈簧的所接受的總虛變形功等于所有彈簧變形功之和,連續(xù)分布彈簧接受總虛功為：

(6)

外力所做的總虛功為：

(7)

由δW外=δW變,將式(6),(7)帶入可得：

(8)

式(6)、式(7)和式(8)中：δδB為虛位移,S和p(x)分別表示內(nèi)部單元作用在邊界單元節(jié)點(diǎn)上的力和作用在邊界單元內(nèi)部的分布力。

由單元剛度方程得：

(9)

然后將式(4)和式(5)帶入式(9)中,得到二維一致粘彈性人工邊界單元的單元剛度矩陣：

(10)

將上面推導(dǎo)過程中的彈簧力用阻尼力替代,使用和上面推導(dǎo)完全相同的方法即可以得到二維一致粘彈性人工邊界單元的阻尼矩陣：

(11)

2 等效二維一致粘彈性邊界單元

2.1 等效二維一致粘彈性邊界單元剛度矩陣

考慮到實際應(yīng)用的方便,尋求一種常用的實體單元等效代替連續(xù)分布的彈簧阻尼元件,為實體單元賦予相應(yīng)的阻尼和彈性模量,實現(xiàn)和連續(xù)分布的彈簧阻尼元件相同的作用,因為這種邊界單元與內(nèi)部單元有相同的形函數(shù),所以認(rèn)為是一致粘彈性邊界單元。采用二維有限元最常用的四節(jié)點(diǎn)矩形單元,來實現(xiàn)這一想法。如圖2所示為四節(jié)點(diǎn)矩形人工邊界單元。

圖1 二維一致粘彈性邊界單元示意圖Fig. 1 Schematic diagram of two-dimensional uniform viscoelastic boundary element

圖2 矩形邊界單元Fig. 2 Rectangular boundary element

令：

(12)

根據(jù)有限單元法基本原理[10],推導(dǎo)平面應(yīng)變時等效剛度矩陣如下：

(13)

為了將四節(jié)點(diǎn)矩形單元和一致粘彈性邊界單元聯(lián)系起來,當(dāng)b趨于0時,可以將式(13)寫成如下形式：

(14)

上飄“～”表示等效參數(shù)。其中：E表示彈性模量,μ表示泊松比。對比式(14)和式(10),兩矩陣完全相等,應(yīng)該有下式成立：

(15)

由式(15)、式(1)和式(2)可推得邊界單元的等效彈性模量和等效剪切模量為：

(16)

2.2 等效二維一致粘彈性邊界單元阻尼矩陣

假設(shè)等效二維一致粘彈性邊界單元阻尼矩陣和剛度矩陣成正比,即：

[Ce]B=β[Ke]B

(17)

其中：

(18)

將式(1)、式(2)、式(10)、式(11)和式(18)代入式(17),可得：

(19)

2.3 邊界驗證算例

由于邊界單元無法分開施加切向阻尼和法向阻尼,本算例中切向和法向阻尼系數(shù)采用表2中幾種不同的組合,進(jìn)行計算比較,并和遠(yuǎn)置邊界、粘彈性邊界和粘性邊界計算結(jié)果進(jìn)行對比。計算工況見表2。

表2 工況表Table 2 Table of working condition

圖3 驗證邊界模型及荷載曲線Fig. 3 Calculation model of verifying boundary and Loading curve

圖4 不同邊界計算結(jié)果對比Fig. 4 Comparison of calculation results of different boundaries

(20)

3 基于子結(jié)構(gòu)的地震動輸入方法

局部人工邊界主要特征是時空解耦,即一個邊界節(jié)點(diǎn)在某一時刻的運(yùn)動僅與其相鄰節(jié)點(diǎn)在相鄰時刻的運(yùn)動有關(guān),而與其他節(jié)點(diǎn)和其他時刻無關(guān)?；谶@個特征,劉晶波等[2]提出了一種適用于各類應(yīng)力型局部人工邊界的地震動輸入方法,通過在包含人工邊界的子結(jié)構(gòu)模型中輸入節(jié)點(diǎn)位移或加速度,經(jīng)過局部子結(jié)構(gòu)動力分析即可獲得相應(yīng)節(jié)點(diǎn)的等效荷載。這種方法可以避免繁瑣復(fù)雜的邊界節(jié)點(diǎn)應(yīng)力的計算。具體實現(xiàn)步驟見圖5。

圖5 子結(jié)構(gòu)法施加地震荷載流程圖Fig. 5 Flow chart of applying seismic load by substructure method

4 不同人工邊界和地震動輸入方法組合的對比

4.1 自由場脈沖荷載輸入

為了驗證人工邊界單元和地震動輸入方法的模擬效果,數(shù)值試驗如下。計算輸入為底部垂直入射的SV位移脈沖,如圖6(b)所示,脈沖持時0.2 s,幅值為0.1 m,計算時長1 s。有限元模型如圖6(a)所示,模型尺寸為40 m×30 m,網(wǎng)格尺寸為1 m×1 m,介質(zhì)密度為2 500 kg/m3,泊松比為0.25,剪切波速為200 m/s,楊氏模量為2.5×108Pa。計算工況見表3。

表3 不同的邊界條件地震動輸入方法Table 3 Ground motion input methods of different boundary condition

圖6 自由場計算模型及輸入位移曲線Fig. 6 Calculation model of free field and Displacement curve

計算結(jié)果如圖7所示。工況3.1與工況3.2計算所得結(jié)果和理論值的擬合度最高。工況3.3同樣能夠得到較為精確的結(jié)果,但是總體精度比工況3.1與3.2得到的結(jié)果精度低,這是因為粘彈性邊界單元在推導(dǎo)時假設(shè)了單元寬度較小并且阻尼系數(shù)使用切向和法向阻尼系數(shù)的平均值。雖然經(jīng)學(xué)者研究[9],粘彈性邊界單元具有較好的魯棒性,但是并不能完全吸收外行波,少量的反射波仍會引起一點(diǎn)誤差。工況3.4采用粘性邊界加子結(jié)構(gòu)法,在于說明其他的人工邊界同樣適用子結(jié)構(gòu)方法輸入荷載,由于將輸入時程基線調(diào)零,并且位移脈沖頻率適中,可以很大程度地減小粘性邊界的漂移現(xiàn)象;工況3.5作為近似人工邊界,在模擬自由場時能夠獲得較高的精度。

圖7 自由場不同點(diǎn)水平位移時程圖Fig. 7 Time history of horizontal displacement at different points in free field

4.2 輸入幅值與邊界單元尺寸比值的影響

由于粘彈性邊界單元實質(zhì)上和普通的有限單元別無二致,當(dāng)應(yīng)變較大時,材料的非線性往往不能忽略,并且在進(jìn)行動力分析時,考慮計算區(qū)域內(nèi)單元材料的非線性經(jīng)常是必要的,在同一分析步中計算的邊界單元也必須考慮過大的應(yīng)變引起的非線性可能導(dǎo)致的計算誤差?？紤]大應(yīng)變可能引起的非線性,再進(jìn)行有限元模型如圖6(a)所示的算例,僅改變輸入位移脈沖的幅值,模型材料參數(shù)與4.1節(jié)相同。以輸入位移幅值umax和粘彈性邊界單元寬度b的比值ε=umax/b為控制邊界變形的參數(shù),分別計算了ε為0.01、0.02、0.05、0.1、0.2和0.5時BE兩點(diǎn)的水平方向的位移反應(yīng),計算結(jié)果如圖8所示。結(jié)果表明：當(dāng)ε>0.1時,BE兩點(diǎn)的計算結(jié)果出現(xiàn)了較大的偏差,并且與理論結(jié)果相差較大,當(dāng)ε≤0.1時,BE兩點(diǎn)位移符合理論解答。表4給出了計算結(jié)果與理論結(jié)果的差值。誤差使用式(21)表示。

表4 BE兩點(diǎn)計算結(jié)果與理論值的誤差Table 4 The errors of the points BE between numerical results and theoretical values

圖8 粘彈性邊界單元不同ε時模型BE兩點(diǎn)水平位移Fig. 8 Horizontal displacement of the model at points B and E with different ε

(21)

綜上所述,說明粘彈性邊界單元的應(yīng)變對其模擬結(jié)果有一定影響,當(dāng)其應(yīng)變較大時,計算結(jié)果出現(xiàn)明顯偏差。建議人為控制ε≤0.1,可以獲得可信度更高的數(shù)值模擬結(jié)果。

4.3 散射場脈沖荷載輸入

建立如圖9所示的散射場模型,驗證幾種人工邊界及地震動輸入方式的組合用于散射場時的計算精度。模型材料的參數(shù)和輸入的位移脈沖和3.2節(jié)完全相同,模型尺寸如圖9所示,計算工況見表3,計算結(jié)果如圖10所示。從工況3.1、工況3.2和工況3.3計算結(jié)果分析,可以得到和3.2節(jié)基本相同的結(jié)論,說明粘彈性邊界和粘彈性邊界單元都有較好的適用性,應(yīng)用于自由場和散射場時都能得到精確結(jié)果。工況3.1與工況3.2和遠(yuǎn)置邊界解的擬合度最高,與3.2節(jié)自由場算例的不同之處在于：(a)工況3.4出現(xiàn)了明顯的漂移現(xiàn)象,這是因為工況3.4采用的粘性邊界,只考慮了對外行波的吸收,可能出現(xiàn)模型整體最終受力平衡時,最終速度不為0的情況,這將導(dǎo)致模型整體沿著最終速度的方向一直運(yùn)動。(b) 工況3.5出現(xiàn)明顯震蕩,原因在于自由度綁定邊界只有在側(cè)邊界掠入射的情況下才能有較高的模擬精度,由于不能保證全部散射波到達(dá)側(cè)邊界時均與側(cè)邊界法線成90°夾角,所以計算結(jié)果出現(xiàn)較大偏差。

圖9 散射場計算模型Fig. 9 Computational model of scattering field

圖10 散射場不同節(jié)點(diǎn)水平位移時程Fig. 10 Time history of horizontal displacement at different points in scattering field

5 結(jié)論

1)粘彈性邊界概念清晰,地震動輸入子結(jié)構(gòu)法節(jié)點(diǎn)力計算準(zhǔn)確,兩者組合計算精度最高,但是在節(jié)點(diǎn)施加彈簧阻尼元件時操作復(fù)雜,前處理工作量大。

2)粘彈性邊界單元和子結(jié)構(gòu)輸入組合模擬精度與集中粘彈性邊界和子結(jié)構(gòu)輸入組合接近,雖然推導(dǎo)粘彈性邊界單元時進(jìn)行了一些假設(shè),但對實際計算精度并沒有很大影響,且有前處理更加簡單方便的優(yōu)點(diǎn)。

3)使用粘彈性邊界單元時,邊界單元的應(yīng)變對計算結(jié)果有一定影響,建議人為控制ε≤0.1,否則計算結(jié)果會產(chǎn)生較大偏差。

4)粘性邊界在輸入地震動頻率合適且輸入時程基線校正時,在自由場計算中能得到比較可靠的結(jié)果,并且在ABAQUS中容易實現(xiàn)。但是應(yīng)用于散射場中時,由于散射波的作用,會出現(xiàn)難以預(yù)估的漂移現(xiàn)象。

5)側(cè)邊界使用自由度綁定邊界時可以精確模擬自由場垂直入射情形,但是對于散射場計算(工況3.5)精度較低,建議在計算散射場時謹(jǐn)慎使用。