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

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

0 引言

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

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

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

1 一致粘彈性人工邊界

1.1 粘彈性邊界理論

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

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

(1)

(2)

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

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

1.2 二維一致粘彈性邊界

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

為了推導(dǎo)一致粘彈性邊界的剛度矩陣和阻尼矩陣,將內(nèi)部單元與粘彈性邊界相鄰的一條邊視作桿件單元,i和j表示單元的兩個(gè)節(jié)點(diǎn),2L表示桿長(zhǎng)。此時(shí),認(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)部的分布力。

由單元?jiǎng)偠确匠痰茫?/p>

(9)

然后將式(4)和式(5)帶入式(9)中,得到二維一致粘彈性人工邊界單元的單元?jiǎng)偠染仃嚕?/p>

(10)

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

(11)

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

2.1 等效二維一致粘彈性邊界單元?jiǎng)偠染仃?/h3>

考慮到實(shí)際應(yīng)用的方便,尋求一種常用的實(shí)體單元等效代替連續(xù)分布的彈簧阻尼元件,為實(shí)體單元賦予相應(yīng)的阻尼和彈性模量,實(shí)現(xiàn)和連續(xù)分布的彈簧阻尼元件相同的作用,因?yàn)檫@種邊界單元與內(nèi)部單元有相同的形函數(shù),所以認(rèn)為是一致粘彈性邊界單元。采用二維有限元最常用的四節(jié)點(diǎn)矩形單元,來(lái)實(shí)現(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)變時(shí)等效剛度矩陣如下：

(13)

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

(14)

上飄“～”表示等效參數(shù)。其中：E表示彈性模量,μ表示泊松比。對(duì)比式(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 邊界驗(yàn)證算例

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

表2 工況表Table 2 Table of working condition

圖3 驗(yàn)證邊界模型及荷載曲線(xiàn)Fig. 3 Calculation model of verifying boundary and Loading curve

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

(20)

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

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

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

4 不同人工邊界和地震動(dòng)輸入方法組合的對(duì)比

4.1 自由場(chǎng)脈沖荷載輸入

為了驗(yàn)證人工邊界單元和地震動(dòng)輸入方法的模擬效果,數(shù)值試驗(yàn)如下。計(jì)算輸入為底部垂直入射的SV位移脈沖,如圖6(b)所示,脈沖持時(shí)0.2 s,幅值為0.1 m,計(jì)算時(shí)長(zhǎng)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。計(jì)算工況見(jiàn)表3。

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

圖6 自由場(chǎng)計(jì)算模型及輸入位移曲線(xiàn)Fig. 6 Calculation model of free field and Displacement curve

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

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

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

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

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

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

(21)

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

4.3 散射場(chǎng)脈沖荷載輸入

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

圖9 散射場(chǎng)計(jì)算模型Fig. 9 Computational model of scattering field

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

5 結(jié)論

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

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

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

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

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