拱壩地震反應(yīng)的能量輻射效應(yīng)

拱壩地震反應(yīng)的能量輻射效應(yīng)

1在地質(zhì)反應(yīng)分析方面的應(yīng)用由于不能忽視拱壩和無限底板之間的能量交換，因此拱壩-地震波前的地震作用下的基質(zhì)系統(tǒng)是能量釋放系統(tǒng)。對于能量開放系統(tǒng)的動力反應(yīng)分析,人工邊界問題是將難以在計算機上實現(xiàn)的無限模型轉(zhuǎn)換為便于在計算機上實現(xiàn)的有限模型的關(guān)鍵。各種常用的人工邊界處理方法參見張楚漢、Wolf等的評述。文獻將無限域的遠場波動位移解假定為無衰減的平面波位移解與無限域格林函數(shù)解的遠場近似位移解的組合,建立了一種新的黏彈性邊界,算例表明,該方法對于近場波源引起的波散射問題模擬具有不低于二階透射邊界的計算精度。拱壩地震反應(yīng)分析方法研究中,Clough基于有限元法開發(fā)了最早的拱壩地震反應(yīng)分析程序——ADAP,其將地基模擬為帶有固定邊界的無質(zhì)量彈性體,沒有考慮無限地基能量輻射影響。Chopra等和Dominguez等在頻域內(nèi)用有限元-邊界元混合法研究了拱壩地震反應(yīng)。Zhang等也提出了拱壩-地基動力相互作用分析的有限元-邊界元-無限元耦合模型。這些頻域求解方法采用的是全局人工邊界條件模擬無限地基能量輻射效應(yīng),存在不便反映結(jié)構(gòu)和近場地基的非線性特征且計算效率低的問題。為提高計算效率,Chopra等和Zhang等采用了近似的阻抗曲線來計算地基施加給拱壩上的作用力,這無疑會影響計算精度。實際上,拱壩-地基系統(tǒng)求解規(guī)模巨大,且往往包含基巖介質(zhì)的節(jié)理裂隙、斷層和軟弱夾層,壩體中縫界面開合,壩體混凝土材料損傷等非線性因素。為了解決這一問題,杜修力等將局部人工邊界中的一種位移型人工邊界——透射邊界引入了拱壩地震反應(yīng)分析。為抑制透射邊界引起的高頻失穩(wěn)現(xiàn)象,DuXiuli等在文獻的分析模型中引進了與介質(zhì)應(yīng)變成正比的阻尼,在一定條件下實現(xiàn)了數(shù)值模型的穩(wěn)定分析,但這一措施并未在理論上圓滿解決透射邊界引起的高頻失穩(wěn)問題。事實上,在使用有限元與透射邊界的結(jié)合法計算拱壩-地基系統(tǒng)的地震反應(yīng)時,為獲得直觀的穩(wěn)定數(shù)值解,一般需要多次反復(fù)試算,這給工程應(yīng)用帶來了不便,也使得可能存在的誤判引起計算精度損失。為克服位移型人工邊界方法存在的這一固有缺陷,本文提出將文獻建立的應(yīng)力型人工邊界方法——一種新的黏彈性邊界用于拱壩-地基開放系統(tǒng)的地震動力反應(yīng)分析,并針對位于Ⅸ度地震設(shè)防區(qū)、壩高292m的小灣拱壩進行地震動力反應(yīng)分析,并與透射邊界方法的分析結(jié)果進行對比。2有限步長的公式2.1人工邊界面的分區(qū)分解拱壩-地基系統(tǒng)內(nèi)結(jié)點有限元運動方程見文獻,無外源作用的黏彈性人工邊界結(jié)點運動方程見文獻。本文進一步給出外源地震作用下黏彈性人工邊界結(jié)點的運動方程及等效地震荷載的計算方法。人工邊界設(shè)置為“盒子”型截斷形式,地震波以平面波形式自底邊界豎直入射,這種邊界截斷和地震動入射形式是多數(shù)結(jié)構(gòu)-地基地震反應(yīng)分析通常采用的。與所有的人工邊界條件一樣,黏彈性邊界僅對從有限域穿過人工邊界進入無限域的外行波的模擬有效,對于入射地震波作用下局部地質(zhì)、地形條件或結(jié)構(gòu)存在引起的地震波散射問題分析,需要在人工邊界處進行波場分解。根據(jù)不同的人工邊界上的波場特點,各邊界面采用兩種不同的分解方案,即底邊界上將總場分解為邊界入射場和邊界外行場,側(cè)邊界上將總場分解為自由場和散射場,邊界入射場或自由場可由連續(xù)介質(zhì)力學(xué)模型解析計算得到,邊界外行場或散射場由人工邊界條件模擬并由離散模型依據(jù)數(shù)值分析方法獲得。這種不同邊界面分區(qū)分解的方法,在盡量避免計算自由場的前提下較好地解決了側(cè)人工邊界面的平行波模擬問題。將未知的散射場或邊界外行場用總場減去自由場或邊界入射場(用上標R表示)表示,得到含外源作用的人工邊界面l結(jié)點i方向的有限元運動方程為mlu¨li+∑k=1ne∑j=1n(clikj+δlkδijAlCli)u˙kj+∑k=1ne∑j=1n(klikj+δlkδijAlKli)ukj=Fli+fli(1)fli=Al(KliuRli+Cliu˙Rli+σRli)(2)mlu¨li+∑k=1ne∑j=1n(clikj+δlkδijAlCli)u˙kj+∑k=1ne∑j=1n(klikj+δlkδijAlΚli)ukj=Fli+fli(1)fli=Al(ΚliuliR+Cliu˙liR+σliR)(2)式中:Kli、Cli為l結(jié)點的人工邊界參數(shù),Kli等于Knlln(i=n,為邊界面法向)或Kτllτ(i=τ,為邊界面切向),Cli等于Cnlln(i=n,為邊界面法向)或Cτllτ(i=τ,為邊界面切向);Al為人工邊界面上l結(jié)點的影響面積;fli為入射地震波作用下l結(jié)點i方向的等效結(jié)點力。式(2)右邊三項分別為克服彈簧、阻尼器及介質(zhì)所需的抗力。底邊界面由入射場直接計算,側(cè)邊界面由自由場計算。2.2等效地震荷載假定在底邊界豎直入射的平面P波和S波的位移時程分別為uP(t)和uS(t),由波動傳播規(guī)律及波場應(yīng)力狀態(tài)按式(2)解析地計算各人工邊界結(jié)點的等效地震荷載fli(t)?？傻萌斯み吔鏻結(jié)點的計算公式如下。(1)定義2.《t1》的材料參數(shù)f?zlz(t)=Al[KnluP(t)+Cnlu˙P(t)+ρcPu˙P(t)](3){f?xlx(t)=AlλcP[u˙P(t?Δt1)?u˙P(t?Δt2)]f?xlz(t)=Al[KτluP(t?Δt1)+Cτlu˙P(t?Δt1)+KτluP(t?Δt2)+Cτlu˙P(t?Δt2)](4){f?yly(t)=?f+xlx(t)=?f+yly(t)=f?xlx(t)f?ylz(t)=f+xlz(t)=f+ylz(t)=f?xlz(t)(5)flz-z(t)=Al[ΚlnuΡ(t)+Clnu˙Ρ(t)+ρcΡu˙Ρ(t)](3){flx-x(t)=AlλcΡ[u˙Ρ(t-Δt1)-u˙Ρ(t-Δt2)]flz-x(t)=Al[ΚlτuΡ(t-Δt1)+Clτu˙Ρ(t-Δt1)+ΚlτuΡ(t-Δt2)+Clτu˙Ρ(t-Δt2)](4){fly-y(t)=-flx+x(t)=-fly+y(t)=flx-x(t)flz-y(t)=flz+x(t)=flz+y(t)=flz-x(t)(5)(2)工作原理f?zlx(t)=Al[KτluS(t)+Cτlu˙S(t)+ρcSu˙S(t)](6){f?xlx(t)=Al[KnluS(t?Δt3)+Cnlu˙S(t?Δt3)+KnluS(t?Δt4)+Cnlu˙S(t?Δt4)]f?xlz(t)=AlρcS[u˙S(t?Δt3)?u˙S(t?Δt4)](7){f+xlx(t)=f?xlx(t)f+xlz(t)=?f?xlz(t)(8)f?ylx(t)=Al[KτluS(t?Δt3)+Cτlu˙S(t?Δt3)+KτluS(t?Δt4)+Cτlu˙S(t?Δt4)](9)f+ylx(t)=f?ylx(t)(10)flx-z(t)=Al[ΚlτuS(t)+Clτu˙S(t)+ρcSu˙S(t)](6){flx-x(t)=Al[ΚlnuS(t-Δt3)+Clnu˙S(t-Δt3)+ΚlnuS(t-Δt4)+Clnu˙S(t-Δt4)]flz-x(t)=AlρcS[u˙S(t-Δt3)-u˙S(t-Δt4)](7){flx+x(t)=flx-x(t)flz+x(t)=-flz-x(t)(8)flx-y(t)=Al[ΚlτuS(t-Δt3)+Clτu˙S(t-Δt3)+ΚlτuS(t-Δt4)+Clτu˙S(t-Δt4)](9)flx+y(t)=flx-y(t)(10)其中:Δt1=l/CP,Δt2=(2L-l)/CP,Δt3=l/CS,Δt4=(2L-l)/C3(11)式中:ρ、cP、cS、λ分別為介質(zhì)密度、P波波速、S波波速、第一拉梅常數(shù);L為底邊界到地表的距離;l為l結(jié)點到底邊界的距離;Δt1和Δt2分別為l結(jié)點處入射P波和地表反射P波的時間延遲;Δt3和Δt4分別為l結(jié)點處入射S波和地表反射S波的時間延遲;等效地震荷載的下標代表結(jié)點號和分量方向,上標代表結(jié)點所在人工邊界面的外法線方向,與坐標軸方向一致為正,相反為負。等效地震荷載fli的計算流程見圖1。2.3第二,是有限元模型的最高固有頻率采用中心差分結(jié)合單邊差分的顯式差分法求解系統(tǒng)運動方程。顯式時間積分是條件穩(wěn)定的,為保證數(shù)值計算收斂,時間步長應(yīng)小于系統(tǒng)臨界時間步長。臨界時間步長是由數(shù)值模型產(chǎn)生的數(shù)學(xué)概念,與模型的最高固有頻率成反比。已經(jīng)證明,由逐個單元為基礎(chǔ)確定的最高單元頻率總是高于有限元模型的最高頻率,這說明由有限個單元構(gòu)成的系統(tǒng)的最高固有頻率是存在的,即有限封閉系統(tǒng)的穩(wěn)定臨界時間步長存在。本文的拱壩-無限地基開放系統(tǒng)通過截斷無限地基施加彈簧-阻尼元件這一黏彈性邊界被轉(zhuǎn)化為形式上封閉的系統(tǒng),因而本文的拱壩—地基系統(tǒng)有限離散模型是存在臨界時間步長的,是條件穩(wěn)定的。3彈性半空間動反應(yīng)問題文獻的研究表明,黏彈性人工邊界計算波源問題(即模擬人工邊界外行波)的精度多數(shù)情況下高于二階透射人工邊界的精度?，F(xiàn)為驗證本文基于黏彈性人工邊界的地震動輸入方法的模擬精度,分析平面P波、S波分別豎直入射下三維彈性半空間的動反應(yīng)問題。入射波位移時程如圖2。介質(zhì)參數(shù)為密度2630kg/m3,彈性模量32.5GPa,泊松比0.22。有限元分析時,截取2000m×2000m×2000m的正方體有限區(qū)域進行計算,采用滿足有限元精度要求的邊長50m的正方體固體單元離散,時間步長0.006s。該問題的自由地表位移解析解為考慮行波延遲后放大2倍的入射位移時程。圖3和圖4分別為P波入射下地表豎向位移反應(yīng)和S波入射下地表水平位移反應(yīng)。比較圖2和圖3及圖4,可以看到,本文方法具有良好的模擬精度。4小灣水庫地震反應(yīng)分析4.1工程區(qū)基本烈度確定小灣水電站位于我國瀾滄江上游河段,總庫容145億m3,總裝機容量4200MW,主要用于發(fā)電、防洪和灌溉,是八級開發(fā)瀾滄江中下游河段的兩個核心電站之一。壩體幾何尺寸和壩體及地基材料特性采用昆明勘測設(shè)計院提供的技術(shù)資料。壩體為拋物線變厚度雙曲拱壩,壩底高程953m,壩頂高程1245m,最大壩高292m,壩頂弧長935m,壩頂厚12m,壩底厚73m,壩體混凝土動態(tài)彈性模量Ed27.3GPa,密度ρd2400kg/m3,泊松比νd為0.189。壩址河谷相對較寬,呈V型,壩頂處河谷寬720m,兩岸山坡平均坡度為40°～42°。壩址區(qū)基巖主要有黑花岡片麻巖和角閃斜長片麻巖組成,細分為21種材料(見表1)。小灣電站周圍地質(zhì)條件比較復(fù)雜,歷史地震發(fā)生頻繁,經(jīng)中國地震局烈度評定委員會審查,確定工程區(qū)基本烈度為Ⅷ度。小灣拱壩的設(shè)計烈度為Ⅸ度,水平向設(shè)計地震動峰值加速度為0.308g,豎向設(shè)計地震動峰值加速度取水平向的2/3,即0.205g。4.2與文獻顯式有限元對比壩體和基巖部分有限元網(wǎng)格采用文獻中的剖分方式,以三維8結(jié)點六面體單元和6結(jié)點退化單元進行有限元離散,共劃分20107個單元,22878個結(jié)點。黏彈性邊界條件采用彈簧-阻尼單元模擬,每一邊界結(jié)點在3個平移自由度方向分別施加彈簧-阻尼單元,一端與人工邊界結(jié)點相連,另一端固定。系統(tǒng)有限元模型如圖5所示。由小灣模型(圖5)看到,壩址河谷的存在使側(cè)向人工邊界面波場較為復(fù)雜,順河向側(cè)人工邊界面包含平行于人工邊界面的自由場和人工邊界外行的河谷散射場,橫河向側(cè)人工邊界面的自由場和河谷散射場均平行于人工邊界面。本文方法的分析結(jié)果將與文獻的顯式有限元結(jié)合透射人工邊界方法的結(jié)果進行對比,三向輸入地震波為與文獻相同的人工合成地震波(如圖6)。計算三向地震作用的無庫水空庫反應(yīng),求解時間增量為0.0005s。計算中采用瑞利阻尼假定,確定瑞利阻尼系數(shù)的兩階阻尼比仍取為0.05(分別對應(yīng)于f1=1.0Hz和f2=15.0Hz)。由于黏彈性邊界和透射邊界等局部解耦的時域邊界條件對平行于人工邊界面的波傳播模擬是無效的,兩種計算方法中分別做如下處理:(1)應(yīng)用透射人工邊界的方法沿用文獻的處理方法,即側(cè)人工邊界面采用透射邊界修正公式,這對于以邊界平行波為主的橫河向側(cè)邊界面而言可以實現(xiàn)完全模擬,但對于不僅存在邊界平行波而且存在大量河谷散射波的順河向側(cè)人工邊界面似乎不盡合理;(2)本文方法中順河向和橫河向側(cè)邊界面的自由場采用式(4)、式(5)、式(7)～式(10)解析確定,順河向側(cè)邊界面散射場由黏彈性邊界模擬,橫河向側(cè)邊界面散射場由黏彈性邊界修正公式模擬。本文處理方法更為合理。4.3結(jié)果分析4.3.1移、速度、加速度時程壩體的最大位移、速度和加速度反應(yīng)為壩頂拱冠順河向反應(yīng),圖7為上游面壩頂拱冠位移、速度和加速度時程曲線。從圖7可以看到,兩種方法的穩(wěn)定階段的結(jié)果基本符合,但隨著計算時間的延長,有限元結(jié)合透射人工邊界的方法的計算結(jié)果出現(xiàn)了高頻失穩(wěn)現(xiàn)象,這表明施加與剛度成正比的阻尼的方法未能根本解決透射邊界引起的高頻失穩(wěn)問題。4.3.2壩面壩面壩面拱冠本文方法計算的壩體的最大主拉應(yīng)力為3.74MPa,最大主壓應(yīng)力為4.55MPa,同樣出現(xiàn)在壩頂拱冠,圖8為上、下游壩面壩頂拱冠主應(yīng)力時程。從圖8可以看到,對于高拱壩這類抗拉強度較弱的大體積素混凝土結(jié)構(gòu),壩頂拱冠是拱壩抗震的關(guān)鍵部位,在高拱壩的抗