地震波斜入射下水平成層覆蓋層地基波動(dòng)輸入研究

王 飛, 宋志強(qiáng), 劉云賀, 李 闖, 李正貴, 胡安奎, 田 慶

(1. 西華大學(xué) 能源與動(dòng)力工程學(xué)院,成都 610039; 2. 西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室,西安 710048;3. 陜西鐵路工程職業(yè)技術(shù)學(xué)院, 陜西 渭南 714099)

隨著壩工技術(shù)進(jìn)一步提高,傳統(tǒng)水電開發(fā)逐步向覆蓋層普遍發(fā)育的西南河流上轉(zhuǎn)移。西南河流河床覆蓋層深厚,其厚度達(dá)數(shù)十米甚至上百米,其中某擬建水電工程壩址覆蓋層最大厚度甚至超過(guò)500 m[1-2],屬超深厚覆蓋層工程。同時(shí),西南地區(qū)強(qiáng)震頻發(fā),地震作用下覆蓋層土體表現(xiàn)明顯非線性特性[3],覆蓋層上大壩抗震安全面臨嚴(yán)峻挑戰(zhàn)。地震動(dòng)輸入是覆蓋層上大壩地震響應(yīng)分析和安全評(píng)價(jià)的基礎(chǔ),目前覆蓋層地基地震動(dòng)輸入大多采用固定邊界結(jié)合一致地震動(dòng)的方法[4-5],即無(wú)質(zhì)量地基輸入模型[6]。無(wú)質(zhì)量地基模型中大壩引起的外行散射波無(wú)法向遠(yuǎn)域逸散,導(dǎo)致大壩地震響應(yīng)被高估[7-8],同時(shí)無(wú)法考慮土體動(dòng)力非線性造成的截?cái)噙吔缟陷斎氲卣鸷奢d的非一致性,該模型難以合理模擬覆蓋層地基地震動(dòng)輸入。

劉晶波等[9]提出的波動(dòng)輸入方法解決了均質(zhì)彈性半空間介質(zhì)輻射阻尼效應(yīng)問題,該方法包括地基截?cái)噙吔缟鲜┘拥酿椥赃吔绾陀勺杂蓤?chǎng)轉(zhuǎn)化而來(lái)的等效結(jié)點(diǎn)荷載,計(jì)算精度高,被廣泛應(yīng)用于基巖上大壩地震響應(yīng)分析[10-12]。覆蓋層土體動(dòng)力非線性和成層特性,引起不同深度自由場(chǎng)和黏彈性邊界參數(shù)在時(shí)間和空間上呈現(xiàn)非線性變化,極大增加了覆蓋層地基自由場(chǎng)和黏彈性邊界參數(shù)求解難度。楊正權(quán)等[13-14]在覆蓋層地基地震動(dòng)輸入中引入黏彈性人工邊界,采用等效線性化方法反映黏彈性邊界參數(shù)隨人工邊界內(nèi)側(cè)土體動(dòng)剪應(yīng)變變化,提高了黏彈性邊界吸能效果,但沒解決覆蓋層地基人工邊界上輸入荷載的問題。Zou等[15-16]建立了人工邊界參數(shù)隨覆蓋層截?cái)噙吔鐑?nèi)側(cè)土體動(dòng)剪應(yīng)變變化的等效人工邊界單元,利用剪切箱模型[17]數(shù)值求解了地震波垂直入射下覆蓋層場(chǎng)地自由場(chǎng),合理解決了地震波垂直入射下覆蓋層地基地震動(dòng)輸入問題。剪切箱模型適用于相同高程運(yùn)動(dòng)一致的情況,地震波斜入射下覆蓋層相同高程土體運(yùn)動(dòng)具有空間非一致性[18],剪切模型不能應(yīng)用求解地震波斜入射下覆蓋層自由場(chǎng)。

劉晶波等[19]提出了平面波斜入射下彈性水平成層半空間自由場(chǎng)一維化有限元時(shí)域計(jì)算法,該方法無(wú)法根據(jù)平坦基巖或覆蓋層表面地震波記錄反演水平成層半空間自由場(chǎng),同時(shí)也未見該方法在覆蓋層地基自由場(chǎng)計(jì)算方面的應(yīng)用。通常平坦基巖或覆蓋層表面地震動(dòng)容易被獲得,基巖或土層中地震動(dòng)難以得到[20]。大多數(shù)情況下需要根據(jù)地表地震動(dòng)反演覆蓋層自由場(chǎng)。

本文引入勢(shì)函數(shù)理論[21],建立覆蓋層頂層波幅矩陣與任意層波幅矩陣之間的動(dòng)態(tài)傳遞關(guān)系,由基巖面入射波和地表面地震動(dòng)分別正演和反演水平成層覆蓋層地基應(yīng)變分量;采用等效線性化方法[22]反映土體動(dòng)力非線性特性,結(jié)合二維平面應(yīng)變理論獲得水平成層非線性覆蓋層地基自由場(chǎng)解析解,進(jìn)而轉(zhuǎn)換為地震輸入荷載。遠(yuǎn)域覆蓋層采用Wang等研究中建立的等效黏彈性人工邊界單元模擬,邊界單元參數(shù)隨土體動(dòng)剪應(yīng)變動(dòng)態(tài)實(shí)時(shí)變化,建立了適用于地震波斜入射下水平成層覆蓋層地基地震動(dòng)輸入方法。通過(guò)兩個(gè)數(shù)值算例驗(yàn)證了本文方法的合理性。

1 水平成層覆蓋層自由場(chǎng)解析計(jì)算

1.1 彈性水平成層覆蓋層

圖1中彈性水平成層覆蓋層由n層土層組成,第i層Lame常數(shù)、剪切模量和密度分別為λi、Gi和ρi,土層厚度為hi。

圖1 地震波斜入射下成層覆蓋層地基分析模型Fig.1 Analysis model of foundation with layered overburden under seismic wave oblique incidence

土體位移運(yùn)動(dòng)如式(1)所示

u=?φ+?ψ

(1)

式中:φ為P波標(biāo)量勢(shì);ψ為SV波矢量勢(shì)。

基于勢(shì)的波動(dòng)方程為

(2)

式中:?2為拉普拉斯算子;cP和cS分別為P波和SV波波速。

第i層中,滿足方程式(2)的勢(shì)函數(shù)為

(3)

第i層x方向和z方向位移可由勢(shì)函數(shù)表示為[24]

(4)

第i層應(yīng)力表示為

將式(3)分別代入式(4)和式(5),得到

(6)

式中:Qi為第i層系數(shù)矩陣;Ei為第i層幅值矩陣;Ui為第i層指數(shù)矩陣。

相鄰?fù)翆臃謱用鏉M足位移和應(yīng)力連續(xù)條件為

(7)

將式(7)代入式(6),可得到

Ei+1=(Qi+1)-1QidiEi

(8)

式中:

根據(jù)式(8)可建立第n層幅值矩陣與第1層幅值矩陣之間的傳遞關(guān)系

(9)

式中,Kn-1,1為第n層和第1層之間的幅值傳遞矩陣。

由式(10)求解最頂層幅值矩陣E1,根據(jù)式(9)中傳遞關(guān)系依次求解每一層幅值矩陣,通過(guò)式(6)經(jīng)傅里葉逆變換獲得土體時(shí)域應(yīng)變分量,也可以獲得位移和應(yīng)力等物理量。上述過(guò)程為基于地表地震動(dòng)直接反演這個(gè)覆蓋層地基線性自由場(chǎng),避免了先反演覆蓋層底部地震動(dòng),再根據(jù)底部地震動(dòng)正演整個(gè)覆蓋層自由場(chǎng)等步驟。

P波入射

SV波入射

求解頂層幅值矩陣E1,然后按照前述過(guò)程求解每一層土體時(shí)域應(yīng)變分量。

1.2 非線性水平成層覆蓋層

(13)

由二維平面應(yīng)變狀態(tài)理論[29],關(guān)系為

(14)

式中:ε1和ε3分別為大、小主應(yīng)變;εxx,εzz和γxz分別為x向、z向正應(yīng)變和剪應(yīng)變,由位移-應(yīng)變關(guān)系可得

(15)

1.3 非線性水平成層覆蓋層自由場(chǎng)解析計(jì)算流程

步驟1依據(jù)實(shí)際地質(zhì)資料建立水平成層覆蓋層地基有限元模型,賦予不同土層相應(yīng)靜力材料參數(shù),計(jì)算得到覆蓋層土體單元靜力初始圍壓;

步驟2將土體動(dòng)力材料參數(shù)和初始圍壓代式(13)計(jì)算土體最大動(dòng)剪切模量,以此作為初始計(jì)算參數(shù),同時(shí)依據(jù)沿深度方向土體最大動(dòng)剪切模量Gmax分布情況對(duì)土層進(jìn)行細(xì)分,即單元最大動(dòng)剪切模量相近的劃分為一小層;

步驟3對(duì)地表地震動(dòng)或建基面入射波進(jìn)行傅里葉變換(fast Fourier transform,FFT),根據(jù)式(10)或式(11)求解最頂層幅值矩陣E1,通過(guò)傳遞矩陣獲得第i層幅值矩陣Ei,經(jīng)傅里葉逆變換(inverse fast Fourier transform,IFFT)獲得土體時(shí)域應(yīng)變分量,利用式(13)計(jì)算不同深度處單元?dú)w一化動(dòng)剪應(yīng)變;

步驟4將歸一化動(dòng)剪應(yīng)變繼續(xù)代入等價(jià)黏彈性模型中,確定新的等效動(dòng)剪切模量和等效阻尼比,按照步驟3重新求解每一小層土體歸一化剪應(yīng)變,以歸一化剪應(yīng)變?yōu)榈諗恐笜?biāo),若本次計(jì)算值與上一次計(jì)算值的誤差滿足允許值,則迭代收斂,否則繼續(xù)迭代,直至迭代滿足允許值;

步驟5將滿足迭代收斂的等效動(dòng)剪切模量和等效阻尼比作為水平成層覆蓋層自由場(chǎng)最終動(dòng)力計(jì)算參數(shù),按照步驟3獲得地震波斜入射下水平成層覆蓋層自由場(chǎng)非線性時(shí)域解析解。

2 地震波動(dòng)輸入方法

地震波斜入射下非線性水平成層覆蓋層地基波動(dòng)輸入方法包括:自由場(chǎng)計(jì)算、動(dòng)力人工邊界參數(shù)計(jì)算和等效輸入荷載計(jì)算。自由場(chǎng)計(jì)算過(guò)程按照1.3節(jié),動(dòng)力人工邊界參數(shù)計(jì)算方法見2.1節(jié),等效輸入荷載由自由場(chǎng)和邊界參數(shù)綜合而得,計(jì)算方法見2.2節(jié)。

2.1 非線性黏彈性人工邊界單元

Wang等研究中建立了黏彈性人工邊界單元,黏彈性人工邊界單元實(shí)現(xiàn)方法:在已有有限元模型外法向延伸一層有限單元,將外層有限單元外側(cè)固定,通過(guò)給外層有限單元賦予等效材料參數(shù)使其作用等價(jià)于集中黏彈性人工邊界。黏彈性人工邊界單元參數(shù)依附于鄰近內(nèi)層單元,通過(guò)等效線性化方法,黏彈性人工邊界單元參數(shù)隨鄰近內(nèi)層單元?jiǎng)蛹魬?yīng)變非線性變化,提高黏彈性人工邊界單元吸能能力。非線性黏彈性人工邊界單元有較高計(jì)算精度,與遠(yuǎn)置邊界計(jì)算結(jié)果相比,最大誤差在5%左右。

黏彈性人工邊界單元等效剪切模量、等效彈性模量和等效泊松比如式(16)、式(17)和式(18)所示

(16)

(17)

(18)

采用與剛度成正比的阻尼集成黏彈性邊界單元阻尼矩陣,阻尼系數(shù)按式(19)取值

(19)

圖2為依附于內(nèi)層土體單元的黏彈性人工邊界單元,非線性黏彈性人工邊界單元實(shí)現(xiàn)流程為:

圖2 黏彈性邊界單元Fig.2 Viscoelastic artificial boundary element

步驟1每次計(jì)算前,由單元震前圍壓、土體動(dòng)力參數(shù)和上次計(jì)算過(guò)程中最大動(dòng)剪應(yīng)變確定單元?jiǎng)蛹羟心Ａ?

步驟2搜索與黏彈性人工邊界單元相同高程且相鄰的土體單元,根據(jù)式(17)、式(18)和式(19)和相鄰單元?jiǎng)蛹羟心Ａ?確定等效彈性模量、等效泊松比和平均剛度比例阻尼系數(shù),賦給對(duì)應(yīng)黏彈性人工邊界單元;

步驟3進(jìn)行下次計(jì)算,獲得地震過(guò)程中新的最大動(dòng)剪應(yīng)變;

步驟4重復(fù)步驟1～步驟3,直至滿足迭代收斂標(biāo)準(zhǔn),得到最優(yōu)的人工邊界單元等效參數(shù)。

2.2 等效輸入荷載

第2節(jié)解析計(jì)算了地震波斜入射下非線性水平成層覆蓋層地基自由場(chǎng),覆蓋層地基截?cái)噙吔缣幾杂蓤?chǎng)轉(zhuǎn)化為邊界單元上等效輸入荷載,實(shí)現(xiàn)地震動(dòng)波動(dòng)輸入。與等效黏彈性人工邊界單元匹配的等效輸入荷載如式(20)所示[30]

(20)

邊界法向

(21)

邊界切向

(22)

在ABAQUS有限元軟件中,在等效黏彈性邊界單元與內(nèi)層土體單元相交面切向和法向施加均布荷載。

2.3 地震波動(dòng)輸入流程

地震波斜入射非線性覆蓋層地基波動(dòng)輸入流程圖,如圖3所示。

圖3 地震波斜入射非線性覆蓋層地基波動(dòng)輸入流程圖Fig.3 Wave input process of foundation with nonlinear overburden under seismic wave oblique incidence

3 水平成層覆蓋層地基波動(dòng)輸入數(shù)值驗(yàn)證

以圖4中彈性水平雙層覆蓋層為研究對(duì)象,土層厚度均為25 m,水平向跨度為100 m,在ABAQUS有限元軟件中采用四結(jié)點(diǎn)雙線性平面應(yīng)變四邊形單元CPE4對(duì)覆蓋層進(jìn)行離散,單元尺寸為1.0 m,表1為彈性水平成層覆蓋層材料參數(shù)。依據(jù)B點(diǎn)入射波和地表面O點(diǎn)地震動(dòng)分別正演和反演覆蓋層底部和側(cè)邊界自由場(chǎng),B點(diǎn)入射波位移時(shí)程為圖5中位移時(shí)程,地表面O點(diǎn)水平和豎向地震動(dòng)位移時(shí)程分別為圖6中位移時(shí)程,采用等效黏彈性人工邊界單元模擬遠(yuǎn)域覆蓋層輻射阻尼效應(yīng)和,將自由場(chǎng)轉(zhuǎn)化為等效黏彈性人工邊界單元表面上輸入荷載,分析波動(dòng)輸入下地表中心A點(diǎn)位移和中部應(yīng)力沿深度變化,并與解析解對(duì)比。依據(jù)B點(diǎn)入射波正演自由場(chǎng)時(shí),P波入射角為30°;依據(jù)點(diǎn)O地震動(dòng)反演自由場(chǎng)時(shí),頂部土層上行P波入射角為30°,頂層土層內(nèi)下行P波反射角,上行SV波和下行SV波入射角和反射角由Snell定律確定,下行P波反射角為30°,上行和下行SV波入射角和反射角均為16.54°。圖7給出了上行P波、下行P波、上行SV波以及下行SV波引起的水平向和豎向位移時(shí)程。

表1 彈性水平雙層覆蓋層力學(xué)參數(shù)

圖4 彈性水平雙層覆蓋層Fig.4 Elastic horizontal double-layers overburden

圖5 正演自由場(chǎng)時(shí)B點(diǎn)入射波位移時(shí)程Fig.5 Displacement time history of the incident wave at point B under forward modeling the free field

圖6 反演自由場(chǎng)控制點(diǎn)O設(shè)計(jì)地震動(dòng)位移時(shí)程Fig.6 Displacement time histories of design ground motion at control point O under inverting the free field

圖7 上行和下行P波以及上行和下行SV波引起的水平向和豎向位移時(shí)程Fig.7 Horizontal and vertical displacement time histories caused by the incident and reflected P-waves, as well as the incident and reflected SV-waves

圖8為依據(jù)B點(diǎn)入射波正演覆蓋層邊界自由場(chǎng)波動(dòng)輸入下表面中心點(diǎn)A位移時(shí)程和中部應(yīng)力峰值沿深度的分布。圖9為依據(jù)地表面點(diǎn)O地震動(dòng)反演覆蓋層邊界自由場(chǎng)波動(dòng)輸入下表面中心點(diǎn)A位移時(shí)程和中部應(yīng)力峰值沿深度的分布。

圖8 正演波動(dòng)輸入下彈性水平成層覆蓋層自由場(chǎng)Fig.8 Elastic horizontal layered overburden free field under forward wave input

圖9 反演波動(dòng)輸入下彈性水平成層覆蓋層自由場(chǎng)Fig.9 Elastic horizontal layered overburden free field under inversion wave input

圖8和圖9表明,地表面中心點(diǎn)A位移時(shí)程和中部應(yīng)力響應(yīng)均與解析解擬合良好,證明了本文建立的波動(dòng)輸入方法即適用于依據(jù)深部基巖中入射波求解地震波斜入射下彈性水平成層覆蓋層地震響應(yīng),也適用于依據(jù)地表地震動(dòng)求解地震波斜入射下彈性水平成層覆蓋層地震響應(yīng),并且計(jì)算精度較高。

以某擬建水電工程壩址河床覆蓋層為原型。根據(jù)地質(zhì)資料,地表以下300 m深度范圍內(nèi),覆蓋層近似水平成層,這里覆蓋層深度取為240 m。如圖10所示,覆蓋層自上而下分為4個(gè)土層,第①層為砂卵石層,厚20 m,第②層為含礫中砂層,厚70 m,第③層為粉質(zhì)黏土層,厚20 m,第④層為粗砂土,厚130 m。首先對(duì)覆蓋層地基進(jìn)行靜力非線性彈性計(jì)算,土體應(yīng)力-應(yīng)變關(guān)系采用鄧肯-張E-v模型模擬,靜力計(jì)算參數(shù)如表2所示,獲取單元震前初始圍壓。土體動(dòng)力計(jì)算采用等價(jià)黏彈性模型,動(dòng)力計(jì)算參數(shù)如表3所示。依據(jù)震前初始圍壓和動(dòng)力計(jì)算參數(shù)確定土體最大動(dòng)剪切模量,對(duì)覆蓋層土層進(jìn)一步細(xì)分,分成33個(gè)小層。

表2 覆蓋層靜力計(jì)算參數(shù)Tab.2 Static calculation parameters of overburden layers

表3 覆蓋層動(dòng)力計(jì)算參數(shù)Tab.3 Dynamic calculation parameters of overburden layers

圖10 非線性水平成層覆蓋層地基Fig.10 Foundation with horizontal layered overburden

為了與該工程場(chǎng)地設(shè)計(jì)地震動(dòng)峰值加速度0.53g匹配,將圖6中水平向和豎向地震動(dòng)幅值分別按2.65倍和1.71倍調(diào)幅,獲得深部B點(diǎn)入射波(取水平向地震動(dòng))和覆蓋層地表面控制點(diǎn)O地震動(dòng)時(shí)程。依據(jù)B點(diǎn)入射波正演自由場(chǎng)時(shí),P波入射角為10°;依據(jù)點(diǎn)O地震動(dòng)反演自由場(chǎng)時(shí),頂部土層上行P波入射角為10°,其余各波型入射角或反射角由Snell定律確定。

圖11和圖12分別為正演和反演非線性水平成層覆蓋層自由場(chǎng)的波動(dòng)輸入下中部位移峰值和應(yīng)力峰值沿深度變化。圖11和圖12均表明,中部位移峰值和應(yīng)力峰值與解析解擬合良好,計(jì)算誤差在5%以內(nèi),計(jì)算精度較高。

圖11 正演波動(dòng)輸入下非線性水平成層覆蓋層自由場(chǎng)Fig.11 Nonlinear horizontal layered overburden free field under forward wave input

圖12 反演波動(dòng)輸入下非線性水平成層覆蓋層自由場(chǎng)Fig.12 Nonlinear horizontal layered overburden free field under inversion wave input

本文建立的波動(dòng)輸入方法適用于依據(jù)地基深部入射波計(jì)算地震波斜入射下非線性水平成層覆蓋層地震響應(yīng),也適用于依據(jù)地表地震動(dòng)計(jì)算地震波斜入射性非線性成層覆蓋層地震響應(yīng),合理解決了地震波斜入射下水平成層覆蓋層地基地震動(dòng)輸入問題。

4 結(jié) 論

地震作用下覆蓋層土體表現(xiàn)明顯的非線性行為,無(wú)質(zhì)量地基模型和通常用于線彈性巖基的波動(dòng)輸入方法難以應(yīng)用于地震波斜入射下非線性水平成層覆蓋層地基。為合理解決地震波斜入射下非線性水平成層覆蓋層地基地震動(dòng)輸入問題,本文主要工作和結(jié)論如下:

(1)建立了地震波斜入射下非線性水平成層覆蓋層地基自由場(chǎng)解析計(jì)算方法。以水平成層覆蓋層地基為研究對(duì)象,基于勢(shì)函數(shù)理論,建立了覆蓋層頂層波幅矩陣與任意層波幅矩陣之間的動(dòng)態(tài)傳遞關(guān)系,得到任意層幅值矩陣,通過(guò)傅里葉逆變換求解了彈性水平成層覆蓋層地基應(yīng)變時(shí)域分量;基于二維平面應(yīng)變狀態(tài)理論和等效線性化方法,迭代計(jì)算了地震波斜入射下非線性水平成層覆蓋層地基自由場(chǎng)時(shí)域解析解;該方法不僅能夠依據(jù)覆蓋層深部斜入射波正演覆蓋層地基自由場(chǎng),還能夠依據(jù)地表面地震動(dòng)反演覆蓋層地基自由場(chǎng),具有避免再根據(jù)地基底部地震動(dòng)正演自由場(chǎng)等過(guò)程的優(yōu)點(diǎn)。

(2)建立了適用于地震波斜入射下非線性水平成層覆蓋層地基地震波動(dòng)輸入方法。利用作者前期發(fā)展的邊界參數(shù)隨邊界內(nèi)層土體單元?jiǎng)蛹魬?yīng)變動(dòng)態(tài)實(shí)時(shí)變化的非線性黏彈性人工邊界單元(發(fā)揮人工邊界最優(yōu)吸能效果),結(jié)合由自由場(chǎng)轉(zhuǎn)換而來(lái)的等效輸入荷載,實(shí)現(xiàn)了地震波斜入射下非線性水平成層覆蓋層地基地震波動(dòng)輸入。研究了覆蓋層深部斜入射波和地表地震動(dòng)已知的兩種情況下彈性水平雙層和非線性水平多層覆蓋層地基地震響應(yīng),計(jì)算結(jié)果均表明,覆蓋層地基位移和應(yīng)力均與解析解擬合良好,驗(yàn)證了本文建立的波動(dòng)輸入方法的高精度性,為地震波斜入射下水平成層覆蓋層地基上土石壩地震響應(yīng)分析提供基礎(chǔ)。