基于諧響應(yīng)的三維場(chǎng)地地震動(dòng)頻域研究

劉爽,閆少敏,裴旺

天津市交通科學(xué)研究院,天津 300300

0 引言

局部場(chǎng)地對(duì)地震動(dòng)存在局部放大或減弱效應(yīng),地形地震動(dòng)響應(yīng)分析需地震學(xué)、工程學(xué)和土木工程學(xué)等基礎(chǔ)科學(xué)輔助[1]。研究局部場(chǎng)地主要有實(shí)測(cè)法、有限元法、邊界元法等方法,其中有限元法的適用性和靈活性較好,應(yīng)用廣泛。

采用有限元法分析局部場(chǎng)地時(shí),為消除人工截?cái)噙吔缟系牟▌?dòng)反射效應(yīng),大大降低有限元分析的計(jì)算自由度,在模擬地基過(guò)程中對(duì)局部地基邊界面添加人工邊界。研究地形效應(yīng)可采用的人工邊界有Sommerfeld邊界、黏性邊界、疊加邊界、Clayton-Engquist邊界、Higdon邊界、雙漸近Higdon邊界、黏彈性人工邊界和透射人工邊界等,黏彈性人工邊界地震動(dòng)分析方法應(yīng)用廣泛。盛志強(qiáng)等[2]采用黏彈性人工邊界的有限元法研究河谷坡角和坡高對(duì)河谷地震動(dòng)響應(yīng)的影響規(guī)律。梁建文等[3]采用有限元軟件ANSYS,基于黏彈性人工邊界理論編寫(xiě)時(shí)域SV波的計(jì)算程序,研究半橢球形和半球形凹陷河谷在平面SV波入射下的地震響應(yīng)。章小龍等[4]將顯示有限元法與黏彈性人工邊界結(jié)合,研究SV波底邊垂直入射作用下復(fù)雜三維地形的地震動(dòng)響應(yīng)。趙瑞斌等[5]采用軟件ABAQUS開(kāi)發(fā)黏彈性人工邊界計(jì)算程序,采用虛擬子空間法研究非對(duì)稱(chēng)階梯地形地震動(dòng)響應(yīng)。巨建民等[6]基于黏彈性人工邊界地震動(dòng)輸入方法理論建立海底沉管隧道模型,分析P波入射下沉管隧道的地震動(dòng)響應(yīng)。賀晨蕊等[7]采用有限元分析法,基于ANSYS分析密集建筑群與沉積谷地在地震荷載作用下的相互影響。Jiang等[8]在賀晨蕊等研究的基礎(chǔ)上引入黏彈性人工邊界,研究三維層狀橢球盆地與城市建筑結(jié)構(gòu)群相互作用的頻域與時(shí)域分析。Ramírez-Len等[9]采用黏彈性人工邊界研究深覆土土石壩動(dòng)力響應(yīng)分析。Liu等[10]結(jié)合ABAQUS,基于黏彈性人工邊界和等效節(jié)點(diǎn)力方法研究3D沉積物的非線性地震響應(yīng)。Wang等[11]采用黏彈性人工邊界理論,對(duì)三維蜿蜒沉積谷底進(jìn)行頻域分析。

采用瞬態(tài)動(dòng)力學(xué)分析模塊研究局部場(chǎng)地效應(yīng)時(shí),瞬態(tài)分析模塊在計(jì)算場(chǎng)地地震荷載時(shí)有極大優(yōu)勢(shì),但荷載步過(guò)多,對(duì)單元數(shù)量為百萬(wàn)級(jí)甚至千萬(wàn)級(jí)的三維場(chǎng)地模型的計(jì)算時(shí)間過(guò)長(zhǎng),不利于程序的糾錯(cuò)及尋找一般性規(guī)律。文獻(xiàn)[12-15]認(rèn)為改進(jìn)黏彈性人工邊界的方法后,可提高計(jì)算精度,未解決計(jì)算效率低的問(wèn)題。基于頻域分析理論[16],線性結(jié)構(gòu)體系中采用疊加理論,工程波動(dòng)領(lǐng)域中時(shí)域結(jié)果與頻域結(jié)果完全等價(jià),基于快速傅里葉變換法,求解頻域結(jié)果較容易[17]。時(shí)域中地震波波形復(fù)雜,荷載隨時(shí)間變化劇烈;在頻域分析過(guò)程中,所施加的荷載是隨時(shí)間變化的簡(jiǎn)諧荷載。諧響應(yīng)分析可研究某一確定頻率簡(jiǎn)諧荷載作用下的結(jié)構(gòu)響應(yīng),可采用諧響應(yīng)分析代替瞬態(tài)動(dòng)力學(xué)計(jì)算分析線性體系結(jié)構(gòu)。在諧響應(yīng)分析過(guò)程中,計(jì)算過(guò)程只有1步,對(duì)隨時(shí)間變化的簡(jiǎn)諧荷載,采用諧響應(yīng)分析求解模塊的計(jì)算時(shí)間遠(yuǎn)少于采用瞬態(tài)分析求解模塊。

本文基于工程波動(dòng)理論,推導(dǎo)平面P波、SV波、SH波垂直入射時(shí)的等效節(jié)點(diǎn)力計(jì)算公式,采用頻域分析相關(guān)理論與快速傅里葉算法提高頻域計(jì)算的計(jì)算效率,將時(shí)域等效節(jié)點(diǎn)力公式轉(zhuǎn)換為頻域計(jì)算公式,通過(guò)自由場(chǎng)地、沉積場(chǎng)地、凹陷場(chǎng)地驗(yàn)證頻域研究方法的準(zhǔn)確性,以期提高三維場(chǎng)地地震動(dòng)響應(yīng)的頻域計(jì)算效率。

1 黏彈性人工邊界與地震動(dòng)輸入方法

圖1 三維黏彈性人工邊界示意圖

黏彈性人工邊界是在人工節(jié)段邊界節(jié)點(diǎn)上施加彈簧器與阻尼器的模擬計(jì)算裝置[18],阻尼器吸收波傳播到邊界處的能量,模擬波在無(wú)限域中的傳播;彈簧器吸收人工邊界的反射波,實(shí)現(xiàn)波動(dòng)在邊界上的透射過(guò)程。

三維黏彈性人工邊界構(gòu)造如圖1所示,剛度和阻尼系數(shù)的計(jì)算公式為:

KBN=αnGAb/R,CBN=ρcpAb,

(1)

KBT=αtGAb/R,CBT=ρcsAb,

(2)

式中:KBN為人工邊界法向彈簧剛度,KBT為人工邊界切向彈簧剛度,CBN為人工邊界法向阻尼系數(shù),CBT為人工邊界切向阻尼系數(shù),αn、αt分別為法向、切向彈簧的修正系數(shù),G為介質(zhì)的剪切模量,R為散射波源到人工邊界結(jié)點(diǎn)的距離,Ab為邊界結(jié)點(diǎn)的影響面積,ρ為介質(zhì)密度,cp為介質(zhì)的P波波速,cs為介質(zhì)的S波波速。

采用黏彈性人工邊界分析局部場(chǎng)地效應(yīng)時(shí),計(jì)算各節(jié)點(diǎn)的等效節(jié)點(diǎn)力,實(shí)現(xiàn)地震動(dòng)輸入。通過(guò)理論推導(dǎo)二維模型的等效節(jié)點(diǎn)力公式[19]。若人工邊界結(jié)點(diǎn)處有自由場(chǎng)位移向量

(3)

式中u、v、w分別為自由場(chǎng)空間沿x、y、z軸方向的位移。

自由場(chǎng)速度向量為

(4)

自由場(chǎng)應(yīng)力張量為

(5)

式中:σxx、σyy、σzz分別為自由場(chǎng)沿x、y、z軸方向的正應(yīng)力,τxy、τyx為自由場(chǎng)沿xOy平面剪應(yīng)力,τxz、τzx為自由場(chǎng)沿xOz平面剪應(yīng)力,τyz、τzy為自由場(chǎng)沿zOy平面剪應(yīng)力。

基于波動(dòng)理論,人工邊界結(jié)點(diǎn)上的等效結(jié)點(diǎn)力

(6)

式中:Kb為人工邊界的彈簧剛度,Cb為阻尼系數(shù),n為邊界外法線方向余弦向量。

實(shí)際計(jì)算時(shí)通常已知地表加速度,采用傅里葉計(jì)算法求解加速度對(duì)應(yīng)的速度和位移,基于一維波理論確定所有自由場(chǎng)的速度和位移。

地震波從模型底面垂直入射自由場(chǎng)的應(yīng)變公式為:

(7)

式中:εxx、εyy、εzz分別為自由場(chǎng)沿x、y、z軸的應(yīng)變,εyz、εzx、εxy分別為自由場(chǎng)沿yOz、xOz、xOy的平面應(yīng)變。

將式(4)代入線彈性材料的應(yīng)力應(yīng)變關(guān)系,可得自由場(chǎng)應(yīng)力

(8)

式中λ為材料梅拉常數(shù)。

由式(1)～(8)可計(jì)算黏彈性人工邊界節(jié)點(diǎn)處等效節(jié)點(diǎn)力。因篇幅原因,考慮模型對(duì)稱(chēng)性,僅推導(dǎo)底邊與其他2組側(cè)面等效節(jié)點(diǎn)力公式,其他側(cè)面計(jì)算方法類(lèi)似。

2 時(shí)域頻域等效節(jié)點(diǎn)力公式

2.1 P波時(shí)域等效節(jié)點(diǎn)力公式

平面P波從底面垂直入射時(shí),在底面邊界上有u=0,v=0,w=w(t),根據(jù)一維波動(dòng)理論有:

w=w0(t-hd/cp)+w0[t-(2H-hd)/cp],

(9)

(10)

(11)

式中:H為地基底面至地基表面距離,hd為人工邊界各節(jié)點(diǎn)到地基底面距離,t為任意時(shí)刻。地基底面邊界上hd=0,n=(0 0 -1)T。

將式(9)～(11)代入式(8),可知底面邊界上自由場(chǎng)應(yīng)力

(12)

將式(9)～(12)帶入式(6)并簡(jiǎn)化后可得邊界荷載為:

式中:等效節(jié)點(diǎn)力的上標(biāo)為人工邊界節(jié)點(diǎn)所在面的外法線方向,與坐標(biāo)軸方向一致時(shí)為正,與坐標(biāo)軸方向相反時(shí)為負(fù);等效節(jié)點(diǎn)力的下標(biāo)為其在x、y、z軸3個(gè)方向上的分量。

同理,可求得其他邊界的等效荷載為:

2.2 SV波時(shí)域等效節(jié)點(diǎn)力公式

平面SV波從底面垂直入射時(shí),在底面邊界上有:u=u0(t),v=0,ω=0,根據(jù)一維波動(dòng)理論可知自由場(chǎng)在任意位置h任意時(shí)刻t滿(mǎn)足

(13)

式中:cs為土層剪切波速,h為黏彈性人工邊界節(jié)點(diǎn)到地基底邊界的距離。

對(duì)式(13)求導(dǎo)可得

(14)

底面邊界上h=0,n=(0 0 -1)T,帶入式(14)可得底面邊界自由場(chǎng)應(yīng)力為:

(15)

將式(15)代入求解,得底邊節(jié)點(diǎn)等效節(jié)點(diǎn)力

可得

同理,可求得其他表面邊界上的等效荷載:

2.3 SH波時(shí)域等效節(jié)點(diǎn)力公式

平面SH波從基巖底面垂直入射時(shí),在底面邊界上有u=0,v=0,v=v0(t),w=0,據(jù)一維波動(dòng)理論有:

底面邊界上h=0,n=(0 0 -1)T,可得底部表面邊界自由場(chǎng)應(yīng)力:

(16)

由式(16)可得底面人工邊界結(jié)點(diǎn)的等效節(jié)點(diǎn)力為:

可得:

同理,可求得其他邊界的等效節(jié)點(diǎn)力荷載:

2.4 頻域等效節(jié)點(diǎn)力公式

在頻域分析中,假設(shè)入射波為正弦波,依據(jù)波動(dòng)理論與頻域計(jì)算相關(guān)理論變換后[16],入射波、反射波的位移分別為:

u1=ei(ωt+kz′),u2=ei(ωt-kz′),

式中:ω為入射波圓頻率,k為波數(shù),z′為節(jié)點(diǎn)深度。

代入相應(yīng)的時(shí)域公式,計(jì)算可得頻域中P波、SV波、SH波各邊界上的等效荷載。

2.5 頻域公式

P波底面邊界上的等效荷載為:

其他邊界上的等效荷載為:

SV波底面邊界上的等效荷載為:

其他邊界上的等效荷載為:

SH波底面邊界上的等效荷載為:

其他邊界上的等效荷載為:

采用軟件ANSYS編寫(xiě)頻域分析的核心程序部分[19]。以P波為例,采用ANSYS諧響應(yīng)分析方法分析場(chǎng)地效應(yīng)的加載與求解,編制頻域分析程序,SV波、SH波的頻域分析程序與P波類(lèi)似。

3 方法驗(yàn)證

3.1 自由場(chǎng)地驗(yàn)證

自由場(chǎng)地驗(yàn)證計(jì)算模型為300 m×300 m×300 m,自由場(chǎng)中心半徑R=75 m。采用軟件Solid45高階實(shí)體單元模擬土體。為體現(xiàn)黏彈性人工邊界良好的適用性,取較高頻率進(jìn)行自由場(chǎng)地驗(yàn)證。計(jì)算參數(shù)為:無(wú)量綱頻率η=2,剪切波速vs=1 km/s,土層密度ρ=2 100 kg/m3,場(chǎng)地阻尼比ζ=0.03,泊松比為1/3。采用簡(jiǎn)諧SV波從底部垂直入射,入射波幅值為1.0。x、y軸上的地表節(jié)點(diǎn)在x方向的位移均約為2.0 m。計(jì)算結(jié)果如圖2、3所示。

a)x軸(y=0,z=0) b)y軸(x=0,z=0) 圖2 x、 y軸上地表節(jié)點(diǎn)x方向位移

a)x軸 b) y軸 圖3 x、y軸誤差分析

根據(jù)經(jīng)典波動(dòng)理論,自由場(chǎng)地表理論位移為2.0 m。由圖2可知:x軸上地表節(jié)點(diǎn)的最大位移為2.01 m,最小為1.92 m,與理論位移最大誤差為4.2%;y軸上地表節(jié)點(diǎn)的最大位移為2.01 m,最小為1.93 m,與理論位移最大誤差為3.7%。由圖3可知,自由場(chǎng)地中各項(xiàng)地表位移與理論位移的誤差小于4.5%。采用諧響應(yīng)分析施加人工邊界的三維場(chǎng)地模型可行,且模型預(yù)測(cè)位移的準(zhǔn)確度較高。

3.2 沉積場(chǎng)地驗(yàn)證

為驗(yàn)證諧響應(yīng)分析求解沉積場(chǎng)地地震動(dòng)響應(yīng)的適用性,以文獻(xiàn)[20]中的半球沉積盆地理論結(jié)論為例,采用半球形沉積盆地作為驗(yàn)證對(duì)象,建立整體模型,在5個(gè)表面加上黏彈性人工邊界進(jìn)行數(shù)值模擬求解。沉積盆地模型的參數(shù)為:沉積盆地半徑R=200 m,沉積盆地內(nèi)剪切波速vsl=500 m/s,壓縮波速vpl=1 km/s,密度ρ=2 100 kg/m3。半空間場(chǎng)地剪切波速vs2=1 km/s,密度ρ2=2 100 kg/m3,沉積場(chǎng)地、半空間自由場(chǎng)泊松比均為1/3,無(wú)量綱頻率分別為0.50、0.75。平面P波、SV波、SH波沿模型底部垂直入射,入射波幅值為1.0。對(duì)比有限元模擬法的計(jì)算結(jié)果與半球沉積盆地理論結(jié)果,如圖4所示。ux、uy、uz分別為x、y、z方向位移。

a)P波 b) SV波 c) SH波 圖4 不同波型入射沉積場(chǎng)地地表位移

由圖4可知:有限元模擬法的計(jì)算結(jié)果與半球沉積盆地理論結(jié)果誤差小于3.0%,準(zhǔn)確度較高,說(shuō)明諧響應(yīng)分析適用于求解沉積盆地地震動(dòng)響應(yīng)問(wèn)題;沉積盆地對(duì)地震波的放大作用更明顯,越靠近盆地中心,放大效果越顯著,沉積盆地對(duì)入射波的聚焦效應(yīng)十分顯著[21-22]。

3.3 凹陷場(chǎng)地驗(yàn)證

為驗(yàn)證諧響應(yīng)分析求解凹陷場(chǎng)地地震動(dòng)響應(yīng)頻域結(jié)果的適用性,以文獻(xiàn)[20]中的半球凹陷場(chǎng)地理論結(jié)論為例,以半球形凹陷場(chǎng)地為驗(yàn)證對(duì)象,凹陷模型的參數(shù)為:凹陷半徑R=200 m,半空間凹陷場(chǎng)地剪切波速vs=1 km/s,密度ρ=2 100 kg/m3,泊松比為1/3,無(wú)量綱頻率η=1.0。平面P波、SV波、SH波沿模型底部平面垂直入射,幅值為1.0。對(duì)比諧響應(yīng)分析計(jì)算結(jié)果與半球凹陷場(chǎng)地理論結(jié)果,結(jié)果如圖5所示。

a)P波 b) SV波 c) SH波 圖5 不同波型入射凹陷場(chǎng)地地表位移

由圖5可知:諧響應(yīng)分析計(jì)算結(jié)果與半球凹陷場(chǎng)地理論結(jié)果誤差小于3.2%,準(zhǔn)確度較高,說(shuō)明諧響應(yīng)分析法適用于求解凹陷地地震動(dòng)響應(yīng)問(wèn)題;平面S波凹陷部分地表的位移低于周?chē)杂蓤?chǎng),平面P波則相反,表明半球形凹陷地形放大或縮小了地表的響應(yīng)。

4 結(jié)束語(yǔ)

基于三維黏彈性人工邊界的局部場(chǎng)地地震動(dòng)輸入方法,根據(jù)波動(dòng)理論,分別推導(dǎo)平面P波、SV波、SH波垂直入射時(shí)的等效節(jié)點(diǎn)力計(jì)算公式。為顯著提高有限元法的計(jì)算效率,采用快速傅里葉變換法將平面SV波、P波、SH波時(shí)域中的等效節(jié)點(diǎn)力荷載轉(zhuǎn)換為頻域中的等效節(jié)點(diǎn)荷載,推導(dǎo)平面SV波、P波、SH波垂直入射時(shí)頻域分析中的等效節(jié)點(diǎn)力計(jì)算公式,并編制ANSYS諧響應(yīng)分析的核心命令流程序。為驗(yàn)證頻域研究方法的準(zhǔn)確性,與自由場(chǎng)地三維沉積場(chǎng)地和三維凹陷場(chǎng)地的結(jié)果進(jìn)行驗(yàn)證對(duì)比。結(jié)果表明,采用頻域等效節(jié)點(diǎn)力公式求解頻域結(jié)果的準(zhǔn)確度較高。

本文僅對(duì)三維垂直入射的正弦波動(dòng)響應(yīng)分析給出頻域計(jì)算公式與相應(yīng)程序,該研究方法與程序適用于二維場(chǎng)地分析。對(duì)其他彈性波(如Rayleigh波、Ricker波等)需進(jìn)一步依據(jù)相關(guān)理論轉(zhuǎn)換,同時(shí)斜波入射的頻域分析仍需進(jìn)一步研究,提高場(chǎng)地效應(yīng)頻域分析的計(jì)算效率。