地震動波動輸入方法在高土石壩動力分析中的應(yīng)用研究

魏匡民 陳生水 李國英 米占寬 沈 婷

(1.南京水利科學(xué)研究院巖土工程研究所，南京 210024；2.水利部土石壩破壞機理與防控技術(shù)重點實驗室，南京 210029)

地震動輸入是大壩地震安全評價的首要前提，也是大壩工程抗震分析中的一個薄弱環(huán)節(jié)[1].目前土石壩有限元動力分析中，大多采用一致性均勻輸入，即將土石壩動力反應(yīng)作為一個能量封閉的振動問題，計算中不計結(jié)構(gòu)與地基的相互作用，將地震慣性力作用在壩體建基面上，則動力計算結(jié)果中加速度、速度、位移均為相對于建基面運動的相對值.一致性地震動輸入方法不能反映河谷不同部位地震動差異和地震動傳播的“行波效應(yīng)”，該方法僅適用于尺寸小，剛度小，輸入頻率低的建筑物[2].對于200～300m 級特高土石壩結(jié)構(gòu)，其尺寸、跨度與質(zhì)量均非常大，在地震分析中應(yīng)考慮壩體與地基的相互作用，這種相互作用既包括半無限地基對壩體動力特性的影響又包括壩體對動輸入的影響.若不考慮壩體與地基的相互作用，則壩體與地基不存在能量交換，尤其不能反映地震能量向無限域的逸散現(xiàn)象，又稱為“輻射阻尼效應(yīng)”.輻射阻尼效應(yīng)和自由場地震動輸入機制密切相關(guān)，地基-壩體動力系統(tǒng)的地震響應(yīng)包含自由場入射地震波以及由壩體和河谷產(chǎn)生的散射外行波，外行波在向無限山體和地基傳播過程中由于幾何擴散和地基阻尼逐步逸散，而實際數(shù)值模擬中計算范圍僅能取1～2倍壩高，若在數(shù)值模型截斷邊界處采用固定約束或自由端，則一部分原本逸散到無限地基中的散射波會經(jīng)截斷邊界反射回壩體-地基系統(tǒng)內(nèi)，顯著影響壩體結(jié)構(gòu)的地震響應(yīng).

為了在數(shù)值計算中合理模擬地基的輻射阻尼現(xiàn)象，眾多學(xué)者開展了截斷邊界處人工邊界的研究，比較有代表性的有：廖振鵬[3]提出的人工透射邊界，Lysmer[4]提出的粘性人工邊界以及以Deeks[5]、劉晶波[6-9]、杜修力[9-12]等學(xué)者為代表的粘彈性人工邊界.相較于粘性邊界和透射邊界，粘彈性邊界在模擬地基的彈性恢復(fù)能力以及低頻穩(wěn)定性方面更有優(yōu)勢[7]，該方法物理概念清晰，數(shù)值穩(wěn)定性高，具有很強的實用性.近些年，粘彈性人工邊界在大壩地震分析中，尤其是混凝土高壩地震分析中，得到了較為廣泛的應(yīng)用，其中杜修力等采用粘彈性人工邊界對小灣拱壩進行了動力反應(yīng)分析[11]，張伯艷等[13]采用粘彈性人工邊界對烏東德拱壩進行了非線性動力分析.在土石壩地震分析領(lǐng)域，孔憲京等[14-15]將非一致性波動輸入方法引入了土石壩動力分析程序，研究了地震波頻譜特性、壩體高度和地基模量對高土石壩地震反應(yīng)的影響.總體來說，非一致波動輸入方法在土石壩中的研究應(yīng)用仍很少，目前土石壩工程規(guī)模不斷增大，一些擬建的大壩如大石峽砂礫石面板壩、如美心墻堆石壩其壩高已達到250～300m 級，土石壩抗震計算中采用更合理的地震動輸入方法有助于更準(zhǔn)確評估大壩實際的抗震性能.本次研究在三維土石壩靜、動力分析程序TOSS3D 基礎(chǔ)上，引入了基于粘彈性人工邊界的非一致性地震動波動輸入方法，通過與解析結(jié)果對比驗證了程序開發(fā)的正確性.以新疆大石峽面板砂礫石壩為例，比較研究了一致性和非一致性地震動輸入下大壩在加速度反應(yīng)、面板動應(yīng)力、震后殘余變形等方面的差異，評價了不同地震動輸入方法對大壩地震反應(yīng)的影響.

1 粘彈性人工邊界單元

人工粘彈性邊界是在地基邊界處人為設(shè)置一種連續(xù)的應(yīng)力型邊界條件，其基本思想是將外行散射波作為波源問題，從理想介質(zhì)中平面、柱面、球面波的標(biāo)準(zhǔn)方程推求結(jié)構(gòu)邊界面上法向和切向應(yīng)力與位移、速度的關(guān)系式，以獲得動態(tài)阻抗矩陣的局部阻尼和剛度集中系數(shù)，相當(dāng)于設(shè)置阻尼器和彈簧的粘彈性邊界.Deeks[5]提出了軸對稱平面應(yīng)變問題中的時域粘彈性邊界.劉晶波等[16]提出了三維靜動力統(tǒng)一人工邊界的簡潔表達，法向彈簧彈性系數(shù)KN，切向彈簧彈性系數(shù)KT，法向阻尼器阻尼系數(shù)CN，切向阻尼器阻尼系數(shù)CT的表達如式：

式 中 ，G為 地 基 剪 切 模 量 ；E為 地 基 彈 性 模 量 ；ρ為 地基 密 度 ；Cs為S波 波 速 ；Cp為P波 波 速 ；rb為 散 射 源到人工邊界的距離；αN，αT分別為邊界法向修正系數(shù)和切向修正系數(shù)，與波源問題中所采用的波型有關(guān)，從地基底部入射的地震波，若其豎向和水平向分量取以P波和S波傳播的平面波時，αN，αT均取0.5.

粘彈性邊界在有限元方法中通過粘彈性邊界單元實現(xiàn)，文獻[6]提出了二維的一致粘彈性人工邊界單元，并推導(dǎo)了二維單元的剛度矩陣、阻尼矩陣，基于Nastran軟件進行了結(jié)構(gòu)的動力分析.本研究在此基礎(chǔ)上開發(fā)了無厚度的三維的粘彈性人工邊界單元，空間坐標(biāo)系中z'方向為單元法向.

局部坐標(biāo)系x'y'z'中，單元結(jié)點位移{ω}、速度與應(yīng)力{σ}之間的關(guān)系可以寫為

由于該單元是無厚度的 下盤結(jié)點1234坐標(biāo)分別與上盤結(jié)點5，6，7，8 重合，單元形函數(shù)可選擇為：

則局部坐標(biāo)中單元剛度矩陣與阻尼矩陣為

其中，應(yīng)變矩陣B和雅克比行列式|J|表達式為

式(4)和(5)均是在局部坐標(biāo)系建立的，實際計算中應(yīng)轉(zhuǎn)換至整體坐標(biāo)系.

2 人工粘彈性邊界地震動輸入

在數(shù)值計算中，需要首先確定半空間自由波場，計算在所有人工邊界上需要施加的荷載.近些年，研究人員推導(dǎo)了地震P波，S波垂直和斜入射時的人工邊界地震動輸入形式[17]，由于地殼介質(zhì)的密度由地表往下隨地層深度而增大，按物理學(xué)中波在不同介質(zhì)中傳播的折射和反射定律，由地殼深部往地表傳播的地震波，其入射方向?qū)⒅饾u接近于垂直水平地表的豎向[2].所以，工程計算中一般采用垂直入射的地震動輸入方式.即將地表面作為半無限體三維上表面，距離表面較遠的平行于表面的一個假想截面，地震波從該面輸入，設(shè)地震波不隨坐標(biāo)X和Y變化，即將地震波作為一維平面波從底部輸入.如圖1所示，在垂直向輸入P波，沿X軸和Y軸方向分別輸入S波.

圖1 地震動輸入的物理模型

本次研究采用文獻[1]給出的各人工邊界地震荷載表達.

3 非一致地震動輸入程序正確性驗證

本次研究基于南京水利科學(xué)研究院開發(fā)的TOSS3D 軟件平臺，實現(xiàn)了人工粘彈性邊界和地震動的波動輸入，為了驗證本次開發(fā)程序正確性，采用經(jīng)典的自由場地基P波，S波傳播案例驗證[10]，計算模型截取邊長2000m 的正方體，模型及有限元網(wǎng)格如圖2(a)所示，粘彈性邊界單元剖分如圖2(b)所示.地基密度2630kg/m3，彈性模量E=32.5GPa，泊松比=0.22，動力計算時取時間步長為 0.006s.

圖2 自由場有限元模型

地震波從三維模型底部垂直輸入，輸入地震波位移、速度和加速度時程表達式如下

位移時程：

速度時程：

加速度時程：

半無限自由場底部輸入位移過程、速度過程、加速度過程分別為u0(t)，u0(t)，¨u0(t)，則任 意 一 點 任意時刻的位移解析解為：

顯然地表位移、加速度解析解為考慮行波延遲后放大2倍的入射位移速度和加速度時程.圖3為P波入射時模型表面中心點(圖1中A 點)的位移、加速度時程，可以看出，模型頂部中心點動位移和加速度幅值是輸入值的2 倍，解析解與數(shù)值解能夠較好吻合.圖4 為S 波入射時模型表面中心點(圖1 中A點)的位移、加速度時程，可以看出，模型頂部中心點動位移和加速度幅值也均為輸入值的2倍，解析解與數(shù)值解能夠較好吻合.可見，本次研究開發(fā)的軟件具有良好的模擬精度.

圖3 P波輸入時地表點解析解與數(shù)值解對比

圖4 S波輸入時地表結(jié)點解析解與數(shù)值解比較

4 大石峽砂礫石面板壩地震分析研究

4.1 工程概況

大石峽砂礫石面板壩壩高247m，正常蓄水位1700m，為擬建的世界第一高面板壩，場區(qū)地震基本烈度7度，攔河壩設(shè)防烈度為8度，設(shè)計地震100年超越概率2%地震動峰值加速度為387gal.典型剖面材料分區(qū)如圖5所示，有限元模型如圖6所示.單元總數(shù)73167個，結(jié)點總數(shù)71554個.填筑和蓄水過程分51個荷載級模擬.

圖5 大石峽砂礫石面板壩材料分區(qū)

圖6 大石峽面板砂礫石壩三維有限元網(wǎng)格

4.2 計算參數(shù)與輸入地震動

數(shù)值模擬中，筑壩料靜力計算模型采用土石壩工程中廣泛使用的“南水”模型[18]，動力模型采用沈珠江動力模型[19]，計算參數(shù)見表1～2，面板和趾板均為C30混凝土.

本次計算采用100年超越概率2%地震動輸入，地震安評部門提供了3條地震動時程，由于并未直接規(guī)定這3條地震動的輸入方向，所以計算過程中有6種荷載組合，計算結(jié)果表明各荷載組合計算壩體動反應(yīng)規(guī)律基本相同，但量值上有差異，本文選最不利荷載組合為例進行分析，該組合中壩軸向、順河向、豎直向分別輸入的加速度時程如圖7所示，其中豎直向加速度峰值為水平向的2/3.一致性地震動輸入計算中地震動從基巖面輸入，非一致性波動地震動計算中，地震動以面力形式輸入.

表1 大石峽筑壩材料南水模型參數(shù)

表2 沈珠江動力模型參數(shù)

圖7 輸入地震動時程

4.3 計算結(jié)果分析

1)加速度反應(yīng)

圖8和圖9分別給出了地震動一致性和波動輸入下壩體的在壩軸向、順河向、垂直向的加速度放大倍數(shù).一致性輸入條件下壩體在軸向、順河向、豎直向的加速度放大倍數(shù)最大值分別為2.83，2.92，3.03，波動輸入條件下，壩體在軸向、順河向、豎直向的加速度放大倍數(shù)最大值分別為2.11，2.70，2.87，可以看出，一致性輸入方法計算的壩體反應(yīng)要大于波動輸入方法.兩種方法計算的壩體加速度極值位置均位于壩體，但分布位置略有不同，從放大倍數(shù)分布還可以看出，采用波動輸入壩頂?shù)摹氨奚倚?yīng)”較一致性輸入弱，此現(xiàn)象與粘彈性邊界吸收了外行波有關(guān).

圖8 一致性地震動輸入加速度放大倍數(shù)

圖9 地震動波動輸入加速度放大倍數(shù)

2)動位移反應(yīng)

圖10和圖11分別給出了地震動一致性和波動輸入下壩體的在壩軸向、順河向、垂直向動位移極值分布.一致性輸入條件下軸向、順河向、豎直向動位移最大值為22.6cm，40.7cm，14.3cm，波動輸入條件下軸向、順河向、豎直向動位移最大值為42.5cm，40.9cm，31.2cm.由于波動方法計算所得的動位移為壩體的絕對位移，一致性輸入方法計算的動位移是壩體相對于基巖面的相對位移，即基巖面的動位移為零.所以，從量值上波動輸入條件下動位移更大，且分布規(guī)律與一致性地震動輸入有較大差異.

圖10 一致性地震動輸入動位移極值

圖11 地震動波動輸入動位移極值

3)震后殘余變形

一致性地震動輸入和波動輸入，大壩殘余變形分布規(guī)律基本一致，即壩軸向殘余變形表現(xiàn)為兩岸向河床方向的擠壓變形，順河向殘余變形基本表現(xiàn)為指向下游，大壩震陷極值出現(xiàn)在壩頂.篇幅所限，圖12給出了波動輸入下壩體的在壩軸向、順河向、垂直向的殘余變形分布.從計算結(jié)果表明，一致性地震動輸入，壩軸向殘余變形指向右岸和左岸最大值分別為10.1 cm 和13.7cm，順河向位移指向下游最大值為62.7 cm，震陷最大值為92.1cm；波動輸入方法壩軸向殘余變形指向右岸和左岸最大值分別為9.5cm 和9.2 cm，順河向位移指向下游最大值為57.1cm，震陷最大值為76.5cm，可見，波動輸入方法計算的壩體殘余變形較一致性輸入小.

圖12 地震動波動輸入壩體殘余變形

4)防滲體應(yīng)力變形

地震過程中混凝土面板內(nèi)會因為地震位移產(chǎn)生一定的動拉應(yīng)力和動壓應(yīng)力.本次計算結(jié)果表明，一致性和波動輸入方法下，面板動應(yīng)力分布規(guī)律基本一致，即順坡向壓應(yīng)力、拉應(yīng)力最大值出現(xiàn)在面板中部偏向壩頂位置，壩軸向壓、拉應(yīng)力最大值出現(xiàn)在面板靠近兩岸偏向壩頂位置.一致性地震動輸入情況下，面板最大軸向動壓、動拉應(yīng)力分別為5.65MPa、5.26 MPa；順坡向動壓和動拉應(yīng)力最大值分別為8.64 MPa，7.98MPa.地震動波動輸入情況下，面板最大軸向動壓和動拉應(yīng)力分別為5.24MPa、3.94MPa；順坡向動壓、動拉應(yīng)力最大值分別為7.44MPa，6.22 MPa.一致性地震動輸入情況下地震期面板周邊縫三向變位分別為錯動39.8mm、沉陷43.4mm、張開21.6mm.地震動波動輸入情況下震后面板周邊縫三向變位分別為錯動38.6mm、沉陷40.8mm、張開19.9mm.

表3為一致性輸入和波動輸入條件下大壩動力響應(yīng)和殘余變形匯總，可以看出，地震動波動輸入面板軸向壓、拉應(yīng)力分別較一致性輸入減小了約7.3%，25.1%，順河向壓、拉應(yīng)力分別減小了13.8%，22.1%.壩體震陷，波動輸入較一致性輸入減小了約17%，地震期周邊縫三向變位錯動、沉陷、張拉分別減小了1.2mm，2.6mm，1.7mm.可見，一致性輸入偏大地預(yù)測了壩體的動力反應(yīng).

表3 波動輸入與一致性輸入計算結(jié)果匯總

5 結(jié) 論

1)波動輸入方法在模擬地基輻射阻尼和地震波傳播的“行波效應(yīng)”方面較傳統(tǒng)的一致性輸入方法更為合理.

2)將粘彈性邊界和均質(zhì)場地地震動輸入方法嵌入了TOSS3D 有限元分析軟件，通過與解析解比較驗證了開發(fā)軟件的正確性.

3)以擬建的247m 大石峽砂礫石面板壩為例，對比研究了地震動一致性輸入方法和波動輸入方法壩體加速度反應(yīng)、動位移、殘余變形以及防滲體應(yīng)力變形方面的差異，結(jié)果表明，由于考慮了半無限地基的輻射阻尼，波動輸入方法計算所得的壩體動力反應(yīng)和殘余變形整體較一致性輸入小，且分布規(guī)律也有所差異.從工程角度來說，傳統(tǒng)的一致性地震動輸入方法是偏于安全的，但地震動波動輸入方法更有助于客觀評價大壩的抗震能力.

4)目前擬建的土石壩壩高已達到250～300m級，由于其質(zhì)量大、跨度長壩體動力分析中需建立壩體與地基的相互作用的系統(tǒng)，以更合理模擬壩體與地基的能量交換行為.