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

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

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

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

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

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

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

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

1.1 彈性水平成層覆蓋層

圖1中彈性水平成層覆蓋層由n層土層組成,第i層Lame常數、剪切模量和密度分別為λ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)的勢函數為

(3)

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

(4)

第i層應力表示為

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

(6)

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

相鄰土層分層面滿足位移和應力連續條件為

(7)

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

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

(8)

式中:

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

(9)

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

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

P波入射

SV波入射

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

1.2 非線性水平成層覆蓋層

(13)

由二維平面應變狀態理論[29],關系為

(14)

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

(15)

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

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

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

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

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

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

2 地震波動輸入方法

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

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

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

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

(16)

(17)

(18)

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

(19)

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

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

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

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

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

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

2.2 等效輸入荷載

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

(20)

邊界法向

(21)

邊界切向

(22)

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

2.3 地震波動輸入流程

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

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

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

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

表1 彈性水平雙層覆蓋層力學參數

圖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為依據B點入射波正演覆蓋層邊界自由場波動輸入下表面中心點A位移時程和中部應力峰值沿深度的分布。圖9為依據地表面點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位移時程和中部應力響應均與解析解擬合良好,證明了本文建立的波動輸入方法即適用于依據深部基巖中入射波求解地震波斜入射下彈性水平成層覆蓋層地震響應,也適用于依據地表地震動求解地震波斜入射下彈性水平成層覆蓋層地震響應,并且計算精度較高。

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

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

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

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

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

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

圖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

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

4 結 論

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

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

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