地震波斜入射下水平成層覆蓋層地基波動輸入研究

王 飛, 宋志強, 劉云賀, 李 闖, 李正貴, 胡安奎, 田 慶

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

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

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

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

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

1 水平成層覆蓋層自由場解析計算

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

土體位移運動如式(1)所示

u=?φ+?ψ

(1)

式中:φ為P波標量勢;ψ為SV波矢量勢。

基于勢的波動方程為

(2)

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

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

(3)

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

(4)

第i層應力表示為

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

(6)

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

相鄰土層分層面滿足位移和應力連續(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)系依次求解每一層幅值矩陣,通過式(6)經(jīng)傅里葉逆變換獲得土體時域應變分量,也可以獲得位移和應力等物理量。上述過程為基于地表地震動直接反演這個覆蓋層地基線性自由場,避免了先反演覆蓋層底部地震動,再根據(jù)底部地震動正演整個覆蓋層自由場等步驟。

P波入射

SV波入射

求解頂層幅值矩陣E1,然后按照前述過程求解每一層土體時域應變分量。

1.2 非線性水平成層覆蓋層

(13)

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

(14)

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

(15)

1.3 非線性水平成層覆蓋層自由場解析計算流程

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

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

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

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

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

2 地震波動輸入方法

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

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

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

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

(16)

(17)

(18)

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

(19)

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

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

步驟1每次計算前,由單元震前圍壓、土體動力參數(shù)和上次計算過程中最大動剪應變確定單元動剪切模量;

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

步驟3進行下次計算,獲得地震過程中新的最大動剪應變;

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

2.2 等效輸入荷載

第2節(jié)解析計算了地震波斜入射下非線性水平成層覆蓋層地基自由場,覆蓋層地基截斷邊界處自由場轉(zhuǎn)化為邊界單元上等效輸入荷載,實現(xiàn)地震動波動輸入。與等效黏彈性人工邊界單元匹配的等效輸入荷載如式(20)所示[30]

(20)

邊界法向

(21)

邊界切向

(22)

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

2.3 地震波動輸入流程

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

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

3 水平成層覆蓋層地基波動輸入數(shù)值驗證

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

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

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

圖5 正演自由場時B點入射波位移時程Fig.5 Displacement time history of the incident wave at point B under forward modeling the free field

圖6 反演自由場控制點O設計地震動位移時程Fig.6 Displacement time histories of design ground motion at control point O under inverting the free field

圖7 上行和下行P波以及上行和下行SV波引起的水平向和豎向位移時程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點入射波正演覆蓋層邊界自由場波動輸入下表面中心點A位移時程和中部應力峰值沿深度的分布。圖9為依據(jù)地表面點O地震動反演覆蓋層邊界自由場波動輸入下表面中心點A位移時程和中部應力峰值沿深度的分布。

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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