SV波斜入射下Nam Ngum 5拱壩-庫水-地基系統的地震響應分析

劉云輝,陳燈紅,潘子悅,趙藝園,劉云龍,陳非凡

(1. 三峽大學 防災減災湖北省重點實驗室,湖北 宜昌 443002; 2. 三峽大學 土木與建筑學院,湖北 宜昌 443002)

0 引言

由于我國處于兩大地震帶之間的特殊地理位置,地質構造十分復雜且勘探條件有限,在遭遇強度大、震源淺、分布廣的地震時往往引發(fā)十分嚴重的地質災害。特別是許多大型混凝土重力壩坐落于與高烈度區(qū)重合的西部地區(qū),因此,對于大壩抗震安全研究是在水利建設中必不可少的一環(huán)。其中利用數值仿真技術來研究大壩抗震問題是非常有效的模擬手段之一。

對于壩體的抗震安全性評價主要有4個重要的環(huán)節(jié):壩址地震動輸入、高壩地震動反應分析、高壩抗震動力安全評價以及壩體和基巖材料動力特性與破壞機理研究[1]。其中,壩體-地基-庫水系統[2-6]地震反應是國內外高壩抗震研究的主要內容。合理的利用地震動輸入形式并建立相應的輸入方法是大壩抗震分析的重要內容[7]。 JOYNER等[8]考慮了彈性基底中的有限剛度,提出等效節(jié)點荷載的理念,很好地解決黏彈性邊界上地震波的輸入。劉晶波等[9-10]將波動輸入方法從一維推廣到二維并給出了一種結構-地基相互作用中地震波動輸入方法。杜修力等[11-13]在拱壩-地基系統的地震響應分析中,為獲得帶有人工邊界的有限離散動態(tài)模型的非耦合方法,將顯式有限元的波動求解方法與黏彈性邊界相結合,通過設置一組人工邊界來模擬出散射波從內部區(qū)域到無限外部區(qū)域的傳播,從而模擬出無限外部區(qū)域。

對于實際的工程而言,地震波的入射形式往往不是如假定的垂直入射,而是以一定的角度輸入到相應的結構體中。在對于混凝土壩等大體積水工建筑物的地震波斜入射研究中,吳兆營[14]以土石壩為依據探究了壩體在地震波斜入射時的最不利輸入。孫奔博等[15]將隱式有限元的求解方式應用于應力型黏彈性人工邊界中,將斜入射問題擴展到三維分析模型中。杜修力等[16]以小灣拱壩為研究對象探究了三維情況下SV波斜入射對拱壩的地震反應的影響。李明超等[17]在利用黏彈性人工邊界分析斜入射下混凝土重力壩的塑性損傷分析時,發(fā)現地基邊界材料的復雜性會影響?zhàn)椥赃吔缣幍刃Ч?jié)點荷載的確定。黃景琦等[18]在文獻中指出大型結構計算單元與節(jié)點數目過多時,通過相關的編程軟件批量處理等效節(jié)點荷載導入至有限元計算軟件中時,每一個節(jié)點施加的等效節(jié)點應力不易于實現。SU等[19]提出了一種在ABAQUS中結構整體靜動力分析中穩(wěn)定切換邊界條件類型的方案,通過對大壩在靜、動荷載作用下的內部和外部波動問題進行了二維有限元分析,驗證了該方法能穩(wěn)定切換邊界條件,且提高計算精度。通過以上相關學者的對比研究中發(fā)現傳統的黏彈性人工邊界方法[20]應用于大型結構的仿真分析時,由于大量的等效節(jié)點力需要通過三方編程軟件處理,使得波動輸入前處理工作量大,也比較復雜,不易于波動輸入方法的實現及工程應用。且在結構的動力響應分析中,因靜、動力分析時邊界條件不同,若同時考慮會使得在進行有限元模型計算時,由于靜力分析與動力分析人工邊界轉換過程中會產生一定的誤差。尤其是當這種方法應用于地震波斜入射的相關研究問題中,如果利用常規(guī)的分析方法進行數值模擬時會使得計算誤差進一步增大,因此需對斜入射問題進行更加深入的研究。

本文在以上學者對地震波輸入方法研究的基礎上提出一種三維黏彈性靜-動力統一人工邊界對地震波斜入射進行研究。并以Nam Ngum 5重力拱壩為研究對象,建立三維拱壩-庫水-地基系統的斜入射模型,推導了任意角度下三維SV波斜入射時的等效節(jié)點荷載,構建了三維SV波斜入射下三維黏彈性靜-動力統一人工邊界,并通過用戶子程序接口編寫了三維SV波的波動輸入子程序,基于大型軟件ABAQUS完成地震動輸入,探究了三維壩體在SV波斜入射作用下的動力響應。

1 三維SV波斜入射模型構建

1.1 黏彈性人工邊界

本文采用劉晶波等[21-22]構建的黏彈性邊界,即在邊界處通過在有限元模型中邊界節(jié)點處建立一系列等效彈性彈簧和阻尼器來實現。如圖1所示,三維黏彈性人工邊界處切向與法向的彈簧剛度與阻尼系數分別為:

圖1 黏彈性人工邊界示意圖Fig. 1 Schematic diagram of viscous-spring artificial boundary

(1)

(2)

式中:cs和cp分別為SV波與P波波速;ρ為介質密度;G為剪切模量;ξ為拉梅常數;r為波源至人工邊界點的距離;α、β為無量綱參數,分別取0.8、1.1。

1.2 基于黏彈性邊界的等效節(jié)點荷載

在杜修力等[23]對地震動輸入方法的相關研究中可以得出,地震波輸入可以通過在黏彈性人工邊界上施加等效節(jié)點力的形式來實現對波運動的轉化,得到黏彈性邊界節(jié)點的切線與法線方向的等效節(jié)點力可以表示為:

(3)

1.3 入射SV波時等效節(jié)點荷載求解

一個入射區(qū)域受到以任意角α1入射的SV波,域內任意一點(x,y)的各項位移時程記為ulx(t),左側人工邊界高度為H,下底邊界和地表面長度為L,如圖2所示。

圖2 三維斜入射示意圖 圖3 三維SV波斜入射反射規(guī)律示意圖Fig. 2 Schematic diagram of three-dimensional oblique incidence Fig. 3 Schematic diagram of oblique incident reflection law of 3D SV wave

左側邊界處任意一點相對于零時刻入射波的波陣面延遲時間分別為:

(4)

式中:Δt1、Δt2、Δt3分別為入射SV波、反射SV波、反射P波相對于零時刻入射波的波陣面的延遲時間。

當SV波以一定的角度斜入射至立方體截斷區(qū)域時,根據波動理論中相關理論可知會產生相應的波形轉換產生相應的SV波和P波,如圖3所示。其中反射角分別為α1、β1,對應的波幅分別為B1、B2。求得左側邊界處的位移場和速度場為:

位移場:

(5)

速度場:

(6)

式(5)與式(6)求解為局部坐標中的位移場與速度場,為求得左側邊界上的應力,則需要將求解出來的位移場與速度場轉化為全局坐標系[24]。

(7)

將左側邊界求得的位移場式(5)與速度場式(6)中入射SV波、反射SV波、反射P波任意一點零時刻相對于波陣面的延遲時間分別替換為式(8)即可按照同樣的方法求得右側邊界下應力:

(8)

前后邊界面應力求解方式與左側邊界求解類似,如圖2中結構所示前側邊界與后側邊界處求解出來的應力為相反數。前側邊界入射的SV波、反射的SV波和反射的P波相對于入射波零時刻波陣面的延遲時間與求解應力為:

(9)

(10)

底邊界面位移場僅由SV波構成,SV波入射波零時刻波陣面的延遲時間與求解應力分別為:

(11)

(12)

1.4 三維黏彈性靜-動力統一人工邊界轉換方案

在結構動力響應分析中,因靜、動力分析時邊界條件不同。在進行整體分析中靜力分析與動力分析之間會涉及邊界條件切換,相關學者[19]研究表明,在邊界條件切換過程中會產生一定的誤差。為消除邊界條件切換產生的誤差對結果產生的影響,本文利用一種“單元生死”的技術處理來實現對大壩-庫水-地基系統邊界條件的轉換。

具體實現的步驟總結如下:①地應力平衡:在動力分析時,為使重力荷載作用時地基仍保持初始狀態(tài),要進行地應力平衡。因此先對結構進行靜力分析,即添加靜力邊界條件和重力荷載,導出地基的應力場。②求解支反力:將上一步導出的應力場作為初始邊界條件,施加靜力邊界條件、重力荷載等靜力荷載進行求解。邊界結點的受力可分解為靜力約束對地基的作用和地基對靜力邊界作用的疊加,最后導出支座反力。③動力分析:將得到的支反力作為靜力約束的等效替代施加于邊界。施加動力邊界條件(黏彈性邊界)以及靜力(包含重力荷載、水壓力等)、地震荷載等進行動力分析。

2 算例驗證

2.1 SV波斜入射算例驗證

本文構造一個三維彈性巖體模型來探究三維SV波斜入射輸入方法的計算精度,如圖4所示,該巖體計算區(qū)域邊長為200 m的立方體塊,網格尺寸取為5 m的立方體,該模型的上方為半空間自由表面,其余邊界全部設置黏彈性邊界。分析時時間增量步為0.01 s,計算時長2 s。位移時程曲線如圖5所示。巖體介質的力學參數如表1所示。入射的SV波位移時間函數如式(13):

圖4 三維彈性巖體有限元模型Fig. 4 Three-dimensional finite element model of elastic rock mass

圖5 入射SV波位移時程曲線圖Fig. 5 Displacement time history of incident SV wave

表1 介質的力學參數Table 1 Mechanical parameters of the medium

(13)

SV波分別在0°、15°及30°斜入射下地表中心點A的水平位移時程曲線如圖6所示,與解析解的位移時程曲線進行對比發(fā)現:SV波斜入射三維算例計算結果的位移時程曲線與解析解的位移時程曲線吻合良好,且?guī)缀醪淮嬖谟嬎憬Y果在位移歸零時出現漂移的情況,可以說明本文中建立的模型及計算程序具有較好的模擬精度。

圖6 斜入射時A點水平位移時程曲線Fig. 6 Horizontal displacements at point A under oblique incidence

2.2 靜動力統一人工邊界算例驗證

為驗證本文中三維黏彈性靜-動力統一人工邊界的模擬精度,選取50 m×50 m×50 m范圍的有限元模型,網格尺寸取2 m的立方體,如圖7所示,其力學參數為:彈性模量E=2.4×1010Pa,泊松比ν=0.25,質量密度ρ=2450 kg/m3,在頂面x向施加一個切向為q=1×108Pa的均布荷載,有限元模型側面和底面施加黏彈性邊界,入射SV波時程曲線與2.1節(jié)中輸入一致,其位移時間函數如式(13),時間步長取0.01 s。

圖7 巖體模型及參考點Fig. 7 Rock mass model and reference point

傳統人工邊界方案中對于靜、動力分析過程中均采用黏彈性邊界來實現結構的動力響應與本文中黏彈性靜-動力統一人工邊界求得頂面中點A處與底面中點B處位移時程曲線圖,如圖8所示,相比于傳統方案,本文中采取的三維黏彈性靜-動力統一人工邊界轉換方案消除邊界切換過程中產生的誤差影響,減小了與理論值之間的誤差,其精度優(yōu)勢能夠得到突顯。說明本文中提出的計算方案與傳統方案相比具有較高的模擬精度。

圖8 模型頂部及底部中點X向位移時程曲線Fig. 8 X displacement time histories of the top and bottom midpoints

3 Nam Ngum 5拱壩地震響應分析

3.1 工程概況

本節(jié)選取Nam Ngum 5重力拱壩為研究對象。Nam Ngum 5重力拱壩的幾何尺寸以及相應的設計資料采用文獻[25]依據,并參照GB 51247—2018《水工建筑物抗震設計標準》[26]確定,最大壩高99.00 m,壩頂寬度6.00 m,拱冠梁底寬42.00 m,心線弧長234.838 m,弧高比2.42,拱頂中心角度92.794°。以此為依據建立Nam Ngum 5重力拱壩-庫水-地基系統的三維有限元模型,如圖9所示,共204538個節(jié)點,190505個單元。其中,壩體有189921個節(jié)點,6440個單元;庫水14616個節(jié)點,12880個單元;地基189921個節(jié)點,171184個單元。

圖9 Nam Ngum 5重力拱壩系統有限元網格圖Fig. 9 Finite element mesh of Nam Ngum 5 gravity arch dam system

3.2 模型材料參數

壩體混凝土考慮為一種混凝土材料,并假定為線彈性,動彈模取為靜彈模的1.5倍[26],阻尼比取5%;地基假定為均質線彈性介質不考慮材料阻尼。材料參數如表2所示,具體壩體、庫水、地基相應的有限元模型以及具體尺寸如圖10所示。

表2 材料力學參數Table 2 Mechanical parameter of materials參數壩體庫水地基靜態(tài)彈性模量/GPa33.15—20密度/(kg/m3)240010002700泊松比0.16—0.30圖10 三維系統示意圖Fig. 10 Three-dimensional system

3.3 計算荷載與邊界條件

對于拱壩-庫水-地基系統的作用荷載主要分為靜荷載與動荷載兩部分,其中靜荷載為自重、靜水壓力等,動荷載為動水壓力與地震荷載。參照規(guī)范可知該地區(qū)地震基本烈度為VI度,地震動峰值加速度為0.15g,本次計算中Δt取0.01 s,地震波是依據規(guī)范譜生成的隨機地震波如圖11所示。

圖11 SV波加速度、位移、速度時程曲線Fig. 11 Acceleration, displacement and velocity of SV wave

采用三角級數法生成人工地震波,對于給定的功率譜密度函數Sx(ω),按式(14)可以方便的生成以Sx(ω)為功率譜密度函數、均值為0的高斯平穩(wěn)過程a(t)[27]。同時為了更好地描述地面運動的非穩(wěn)定性,采用包絡函數f(t)乘以平穩(wěn)過程a(t),具體如式(15)所示的人工地震波型。

(14)

x(t)=f(t)·a(t)

(15)

在系統模型中,截斷邊界采用黏彈性邊界以模擬無限地基的輻射阻尼效應,采用聲學單元及無反射邊界條件模擬壩-庫動力相互作用[28]。在進行結構整體靜動力分析中,通過一種“單元生死”技術,實現靜態(tài)約束邊界到動態(tài)黏彈性邊界的穩(wěn)定切換,從而實現靜動力統一計算。

4 結果分析

利用合成的三維SV波實現對Nam Ngum 5拱壩-庫水-地基系統的位移、應力、動水壓力的分析,由于SV波入射時存在臨界角[29]且公式為α=arcsin(cs/cp)=32.3°,故本次選取0°、5°、10°、15°、20°、25°來進行工程實例分析。

4.1 位移分析

壩頂點橫河向與順河向位移曲線,如圖12所示。由圖可知,隨著入射角度的增加橫河向位移在逐漸減小,25°較 0°斜入射減小了8.32%,順河向位移在逐漸增加,其中在地震波25°斜入射下順河向位移達到0.15 m,這主要是因為輸入的地震荷載為SV波(橫河向位移)且本文中的地震動輸入僅僅考慮平動分量,使得地震波在地表發(fā)生了反射使得地表結構產生了順河向位移,導致水平位移逐漸減小順河向位移逐漸增大。

圖12 不同角度斜入射下壩頂點橫河向、順河向位移Fig. 12 Cross-river and along-river displacements with different angles of oblique incidence at the dam crest

4.2 應力分析

不同入射角度時壩體最大與最小主應力等值線圖如圖13和圖14所示,由圖可知,上游面主拉應力與下游面主拉應力均成對稱分布,且隨著入射角度的增大壩體主拉、壓主應力逐漸增大,壓應力極值集中出現在壩趾處,拉應力極值集中在壩踵處,規(guī)律性均與實際相符。分別提取壩踵、壩趾2個關鍵點的應力峰值曲線如圖15所示。隨著入射角度的逐漸增大,壩踵點、壩趾點2個關鍵點的應力均呈現增大的趨勢,且在地震波25°斜入射下相應的主拉應力與主壓應力達到了1.47 MPa與3.00 MPa。而入射角度為0°(即垂直入射)與其他入射角度相比,最大或最小主應力是最小的,這說明三維情況下SV波隨著入射角度的增加對壩體動力響應的影響越大。

圖13 SV波斜入射時壩體最大主拉應力等值線圖Fig. 13 Contour map of the maximum principal tensile stress in the dam at oblique incidence of SV waves

圖14 SV波斜入射時壩體最大主壓應力等值線圖Fig. 14 Contour map of the maximum principal compressive stress of dam at oblique incidence of SV wave

圖15 壩體主應力峰值Fig. 15 Principal peak stress of dam body

4.3 動水壓力分析

壩踵點處動水壓力時程曲線如圖16所示,由圖可知:在相同位置隨著入射角度的增加動水壓力也在呈現出逐步增大的線性關系,當入射角度為25°時達到最大,此時最大動水壓力峰值為0.33 MPa,最小動水壓力峰值為-0.29 MPa。入射角度為0°時最小,此時最大動水壓力峰值為0.05 MPa,最小動水壓力峰值為-0.044 MPa。這說明三維情況下SV波隨著入射角度的增加對大壩的動水壓力也會逐漸增大。

圖16 壩踵點動水壓力Fig. 16 Hydrodynamic pressure at the dam heel

5 結論

本文將黏彈性人工邊界與相應地震動輸入方法相結合,推導了任意角度下三維SV波傾斜入射的等效節(jié)點荷載,構建了地震波斜入射的三維黏彈性靜-動力統一人工邊界,并且通過對ABAQUS的二次開發(fā),編寫了三維SV波斜入射的等效節(jié)點荷載子程序;構建了三維SV波斜入射拱壩-庫水-地基系統分析模型,采用所建立的傾斜入射的模型與分析方法,將其應用于Nam Ngum 5重力拱壩的地震響應分析中,實現對模型中壩體在地震波傾斜入射波時的位移、應力、動水壓力的分析,得到以下幾點結論:

1)三維SV波斜入射與垂直入射下相比,重力拱壩壩體關鍵點地震響應規(guī)律存在明顯差異。且SV波斜入射時壩體關鍵點的順和向位移、主應力及動水壓力與入射角度呈正相關,其中25°入射情況下的壩體關鍵點位移、主應力、動水壓力最大。

2)由于輸入的地震荷載為SV波(橫河向位移),地震波在地表發(fā)生了反射使得地表結構產生了順河向位移,導致水平位移逐漸減小順河向位移逐漸增大。其中25°入射與0°入射相比,隨著入射角度的增加橫河向位移呈現逐漸減小的趨勢,順河向位移在逐漸增加,其中在地震波25°斜入射下順河向位移達到最大,為0.15 m。

3)隨著入射角度的增加壩體關鍵點的動水壓力及主應力均在不斷增加,其中25°與0°斜入射相比,在地震波25°斜入射下主拉應力與主壓應力達到最大分別為1.47 MPa與3.00 MPa;動水壓力也是25°入射時達到最大,最大動水壓力峰值為0.33 MPa,最小動水壓力峰值為-0.29 MPa。這表明三維SV波斜入射情況下,隨著入射角度的增加對大壩的主應力與動水壓力均會產生一定的影響。