林皋

（大連理工大學(xué)，遼寧省大連市 116024）

大壩抗震分析與安全評價

林皋

（大連理工大學(xué)，遼寧省大連市 116024）

基于大連理工大學(xué)近期開展的一些研究工作，對大壩抗震分析和安全評價有關(guān)的一些重要問題提出了看法和討論意見。主要包括壩—庫水動力相互作用、壩—地基動力相互作用、大壩抗震分析的精細化數(shù)值模型、大壩的超載潛力與抗震安全性評價。

比例邊界有限元方法；大壩抗震特性；壩—庫水動力相互作用；壩—地基動力相互作用；高精度應(yīng)力分析

0 引言

現(xiàn)代意義上的壩工設(shè)計理論和技術(shù)實踐出現(xiàn)在20世紀以后，迄今只有100年左右的歷史。20世紀30年代美國建設(shè)221m高的胡佛重力拱壩，提出了大體積混凝土建筑的溫控冷卻技術(shù)和拱壩設(shè)計的試載法，對混凝土大壩的建設(shè)起了重要的促進作用。第二次世界大戰(zhàn)以后的恢復(fù)期，20世紀70年代，適應(yīng)各種復(fù)雜地形、地質(zhì)條件的新型式的拱壩、支墩壩和連拱壩，在歐洲、蘇聯(lián)和日本等地興建，工程造價大量節(jié)省，建壩速度顯著加快。這一時期，土石壩的設(shè)計理論和筑壩技術(shù)也有了很大進展。美國奧洛維爾土石壩（高230m）、加拿大的邁卡土石壩（高242m）分別在1968年和1972年建成，刷新了當時200m高度土石壩的記錄。蘇聯(lián)（現(xiàn)塔吉克斯坦）1980年建成努列克壩（高300m），則創(chuàng)造了300m高土石壩的記錄。

改革開放和經(jīng)濟建設(shè)的高速發(fā)展，中國高壩建設(shè)也迎來了燦爛輝煌的年代。進入21世紀前后，依靠自力更生，一大批接近和超過世界頂級高度300m的大壩已在我國建成和正在開工建設(shè)，世界建壩的重心已經(jīng)轉(zhuǎn)向中國。大壩建設(shè)的成就也促進了我國大壩抗震技術(shù)的發(fā)展。通過科技攻關(guān)和國家自然科學(xué)基金重點項目的支持，我國大壩抗震研究水平已進入世界先進行列。在強地震活動區(qū)我國已建成了高度居世界前三位的錦屏一級（高305m）、小灣（高294.5m）和溪洛渡（高285.5m）高拱壩，設(shè)計地震加速度分別達到0.269g、0，313g和0.321g。并在設(shè)計地震加速度高達0.557g的高強地震區(qū)建成了高210m的大崗山高拱壩。我國建設(shè)中的雙江口心墻堆石壩，高314m，設(shè)計地震加速度0.208g，建成后，將成為世界上最高的大壩。在2014年印度尼西亞召開的國際大壩委員會82屆年會期間，中國大壩工程學(xué)會將介紹我國大壩抗震研究成果的單行本分發(fā)各國大壩委員會主席、秘書長時，得到各國委員會高度評價和贊揚。相比起來，世界上已建成運行的位于地震活動區(qū)的300m級的高拱壩只有蘇聯(lián)（現(xiàn)格魯吉亞）的英古里拱壩一座，高271.5m，設(shè)計地震加速度0.23g。在歐洲，只有瑞士1957年建成的莫瓦桑拱壩一座，高250m。在抗震研究開展富有成果的國家中，美國除1936年建成的221m高的胡佛重力拱壩外，高度超過200m左右的大壩為數(shù)較少，日本還沒有建設(shè)200m以上的大壩。

我國在以拱壩為代表的大壩抗震理論和技術(shù)方面取得了不少有意義的成果，可成功地求解大壩—庫水動力相互作用、大壩—無限地基動力相互作用、混凝土材料動態(tài)特性、橫縫動接觸計算模型、大壩非線性地震響應(yīng)、大壩地震損傷與斷裂的數(shù)值模擬、大壩模型動力試驗等復(fù)雜課題，并在大壩抗震響應(yīng)的計算精度與計算效率方面與國外技術(shù)相比毫無遜色，甚至優(yōu)于國外技術(shù)。雖然取得了這樣一些重要進展，但是也要看到，在大壩抗震的計算理論與規(guī)范標準方面，我國仍然沒有脫離美國研究者所形成的一些框架，距離引領(lǐng)世界大壩抗震研究的水平還存在一定差距，有待于我國科技人員的共同努力。

本文將根據(jù)大連理工大學(xué)在大壩抗震研究方面所開展的工作，對大壩抗震研究的關(guān)鍵技術(shù)問題、目前所取得的進展以及需要進一步深入探討的問題談一點認識和看法，供討論和參考。

1 大壩抗震安全的前沿理論和技術(shù)

大壩抗震安全的前沿理論和技術(shù)問題主要體現(xiàn)在以下三個方面：

（1）場地地震風險的評估，據(jù)此提出恰當?shù)拇髩慰拐鹪O(shè)防參數(shù)；

（2）大壩的計算模型和計算方法的改進，以便對大壩地震響應(yīng)做出比較準確而可靠的預(yù)測；

（3）大壩抗震裕度和超載潛力的分析，以便對大壩的抗震設(shè)防的有效性做出較為可靠的評估。

需要指出，這些問題目前并沒有很好地得到解決，特別是第一個和第三個問題，研究工作基本上還處在發(fā)展階段。各個國家的看法和做法還各不相同，也各有特色。值得指出，有關(guān)大壩地震響應(yīng)和震害的實測資料還十分缺乏，經(jīng)受過強震考驗的大壩為數(shù)甚少，只對高度為百米級的混凝土壩和土石壩獲得一些不完整記錄或震害觀測，包括美國的帕柯依瑪（Pacoima）拱壩（高113m）、印度的柯依那（Koyna）重力壩（高103m）、我國的新豐江大頭壩（高105m）、伊朗的瑟菲特·盧德（Sefid Rud）大頭壩（高106m）、我國的紫坪鋪面板堆石壩（高156m）和沙牌拱壩（高132m）等，這些情況也增加了大壩研究的復(fù)雜性和難度。

現(xiàn)有的大壩抗震安全評價主要基于線彈性動力分析。著名壩工專家Chopra在2008年北京召開的第14屆世界地震工程會議[1]和1992年西班牙馬德里召開的第10屆世界地震工程會議的特邀報告[2]中指出，拱壩的抗震分析應(yīng)考慮以下各種因素：水庫和地基的無限性、壩—水相互作用、水庫邊界對波的吸收作用、水的可壓縮性、壩—基相互作用、壩基的地震動不均勻輸入的影響。他認為現(xiàn)有的一些大壩抗震設(shè)計方法比較簡單而不夠準確，沒有全面考慮這些因素、Chopra的觀點值得研究。

根據(jù)上述情況，下文將對大壩抗震分析和安全評價的一些重要問題分別進行探討。

2 壩—庫水動力相互作用

壩—庫水—地基體系的地震響應(yīng)分析（見圖1），需要考慮壩—庫水動力相互作用以及壩—地基動力相互作用的影響。壩—庫水的動力相互作用問題，長期以來受到研究者的重視。自1992年Westergaard提出基于勢流理論的動水壓力計算方法以來，已經(jīng)有大量論文發(fā)表。根據(jù)地震波特性和振動頻率所做的某些基本假定，研究者論述了忽略水的可壓縮性的附加質(zhì)量模型的計算結(jié)果與考慮庫水可壓縮性的計算結(jié)果相近。于是，以上游直立面重力壩剛性壩面動水壓力為基準的附加質(zhì)量計算公式在壩工設(shè)計中得到廣泛的推廣應(yīng)用，并且一直沿用至今。

式中：pw——水深h處的動水壓力強度；

ah——水平地震加速度；

ρw——庫水質(zhì)量密度；

H——壩前庫水深度。

實踐表明，式（1）是非常保守的，所以有的工程在設(shè)計中引入了0.5的折減系數(shù)，這實際上很難反映剛性壩面動水壓力與彈性振動動水壓力之間的差別。我國水工建筑物抗震設(shè)計規(guī)范在考慮壩彈性振動的影響后，對式（1）給出的動水壓力隨水深變化的系數(shù)，依據(jù)十幾座重力壩在若干典型地震波作用下動水壓力的統(tǒng)計資料，作了必要調(diào)整。

圖1 壩—庫水—地基系統(tǒng)的地震響應(yīng)Fig.1 Earthquake response of dam-reservoir-foundation system

20世紀80年代起，Chopra團隊開展了壩—水庫動力相互作用，特別是拱壩動水壓力有限元計算方法的系統(tǒng)研究[3][4]，加深了對動水壓力的理解。Chopra根據(jù)對若干拱壩所作計算分析的結(jié)果指出，不考慮庫水的壓縮性可能顯著低估拱壩動水壓力的大?。ɡ鏜onticello拱壩），也可能顯著高估動水壓力的大?。∕orrow point 拱壩），但Chopra提出的拱壩動水壓力計算方法比較復(fù)雜，工作量大，所以在壩工設(shè)計中難以推廣應(yīng)用。

我們提出了基于比例邊界有限元法（SBFEM）的大壩動水壓力的計算方法[5]，SBFEM是一種高效的半解析半數(shù)值的計算方法，特別適用于具有無限域問題的求解。SBFEM對三維水庫只需進行計算域的表面離散，簡化為二維問題，對于棱柱形水庫（水庫斷面沿軸線不發(fā)生變化），只需在壩面進行離散（見圖2），計算工作量大為節(jié)省。SBFEM可以很方便地考慮庫水壓縮性和水庫邊界對波的吸收作用。

圖2 拱壩動水壓力計算壩面離散網(wǎng)格Fig. 2 Mesh discretization along dam-reservoir interface for hydronamic pressure analysis of arch dam

動水壓力的計算式表示如下：

式中：

符號{ph}、{pc}、{pv}分別代表水平順河向、水平橫河向（或豎向）地震時作用于壩面的動水壓力；代表地震產(chǎn)生的壩面法向加速度，計算ph時，代表順河向地震產(chǎn)生的壩面法向加速度，計算pc時，ün則代表橫河向地震或者豎向地震產(chǎn)生的壩面法向加速度。

[Ф]和[?]代表水庫的振動模態(tài)與振動頻率，滿足方程（6）和壩面的幾何形狀與網(wǎng)絡(luò)離散特性（以哈密頓矩陣[Z]表示）相關(guān)。

計算式（2）～式（8）表明，我們將動水壓力表示為庫水振動模態(tài)和振動頻率的函數(shù)，給出了顯式表達式；Chopra提出的有限元方法將動力水壓力表示為壩體振動模態(tài)的函數(shù)，計算復(fù)雜。相比起來我們的方法計算簡便、精度高，而且收斂性好。

我們計算了美國Morrow Point拱壩壩面動水壓力的頻響函數(shù)，并與Chopra等的計算結(jié)果相比較，結(jié)果如圖3所示，解的相符性很好。圖中曲線表明水庫岸坡和庫底的吸收作用（以反射系數(shù)α表示）對動水壓力的幅值影響較大。

據(jù)此，我們計算了Morrow Point拱壩（大壩假設(shè)為對稱形狀，可以只計算半個壩）在Koyna地震作用下（峰值加速度假設(shè)為0.2g）產(chǎn)生的第一主應(yīng)力幅值沿壩面的分布，計算中考慮了壩—庫水動力相互作用和壩—地基動力相互作用的影響。為了比較，動水壓力采用了不同的計算模型。除所建議的方法SBFEM外，還采用了附加質(zhì)量方法（忽略庫水壓縮性，動水壓力滿足拉普拉斯方程，附加質(zhì)量沿壩面分布為滿陣）以及按經(jīng)典Westergaard公式[按式（1）計算，附加質(zhì)量分布為對角陣]計算的方法，見圖4。由圖可見，按經(jīng)典的Westergaard公式計算，給出的應(yīng)力計算值最為保守，偏大很多。

這只是問題的一個方面。需要特別指出的是，考慮庫水的可壓縮性與庫岸、庫底對波的吸收作用以后，壩面地震動水壓力與壩體地震慣性力之間存在相位差，壩面各點地震動水壓力之間也存在相位差[5]。這就是說，壩面各點地震動水壓力的最大值不在同一瞬間發(fā)生，與壩體慣性力之間也存在時差。從而地震動水壓力對大壩地震響應(yīng)的影響就要緩解得多。這是很值得重視的。

圖3 Morrow Point拱壩動水壓力頻響函數(shù)Fig.3 Frequency response function for hydrodynamic pressure acting on Morrow Point arch dam

圖4 Morrow Point拱壩在Koyna地震作用下（a=0.2g）采用不同的動水壓力計算模型引起的第一主應(yīng)力沿上、下游壩面的分布Fig.4 Distribution of first principle stresses over upstream and downstream face of Morrow point arch dam subjected to Koyna earthquake excitation（0.2g）by applying various models of hydrodynamic pressure

從以上的分析可見，考慮水的可壓縮性及庫岸、庫底對波的吸收作用以后，計算的動水壓力幅值較附加質(zhì)量法及Westergaard經(jīng)典公式的結(jié)果降低很多。引入庫岸庫底反射系數(shù)后，動水壓力幅值將進一步降低。再考慮到動水壓力幅值與壩體慣性力幅值間存在的相位差，以及動水壓力各點幅值間也存在相位差的因素，實際的地震動水壓力將會降低很多。

以上指出對于棱柱形水庫動水壓力的計算，采用SBFEM分析時不必進行水庫的離散。但對幾何形狀不規(guī)則的水庫，則需將水庫劃分為近區(qū)和遠區(qū)兩部分。幾何形狀不規(guī)則的近區(qū)采用SBFEM進行離散（也可以用FEM離散），遠區(qū)可簡化為形狀規(guī)則的水域以模擬水庫無限域的影響，在近區(qū)和遠區(qū)之間設(shè)置波能量傳遞邊界TBC（見圖 5）。

圖5 幾何形狀規(guī)則水庫與幾何形狀不規(guī)則水庫TBC的設(shè)置Fig.5 Geometrically regular and irregular reserroir and placing the TBC

文獻中廣泛采用的傳遞邊界模型有Sommerfeld輻射邊界模型與Sharan的改進模型[6]。我們研究發(fā)現(xiàn)Sommerfeld TBC與Sharan TBC對于反映水庫無限域的輻射特性方面都存在一定的問題，說明如下。

以圖5（a）所示具有無限水域的重力壩動水壓力計算為例。采用有限元FEM方法分析時，水庫域需要離散，在一定距離L處進行截斷，設(shè)置傳遞邊界。采用SBFEM分析時，對于形狀不變的水庫可以不必離散，相應(yīng)于L=0。但對于圖5（b）所示的不規(guī)則水庫，則需在水庫近區(qū)進行離散，在一定距離處截斷設(shè)置TBC。為了對各種計算模型的效果進行對比，假設(shè)重力壩水庫深度H=100m，各種計算模型都在L=100m處設(shè)置TBC。據(jù)此計算出壩面動水壓力的頻響曲線，并與解析解進行對比，如圖6所示。

由圖可見，在水庫第一振動頻率附近，Sommerfeld TBC與SharanTBC的結(jié)果都發(fā)生較大的偏離，只有我們提出的SBFEM TBC與解析解完全吻合。但當TBC設(shè)置在L=1800m處時，Sommerfeld TBC與Sharan TBC才與解析解相吻合，如圖7所示。也就是說采用Sommerfeld TBC或者Sharan TBC時，對水庫必須進行比較大范圍的離散，才能使計算出的壩面動水壓力獲得較為正確的結(jié)果。

2013年10月2-4日國際大壩委員會在奧地利Graz舉辦了第12屆大壩數(shù)值分析的基準研討會，會上一個重要課題就是地震作用下壩—庫水的動力相互作用分析。8個國家的代表提交了11篇論文[7]，共同對220m高的一座拱壩進行了分析?；疽恢碌慕Y(jié)論是采用Westergaard經(jīng)典公式計算動水壓力過于保守，并使計算的拱壩第1、2階基本頻率偏低，根據(jù)各家研究結(jié)果的統(tǒng)計得出的拱壩第1、2階的基頻約為1.54Hz，但采用Westergaard經(jīng)典公式計算動水壓力時，第1、2階基頻可降為1.26Hz和1.32Hz?？紤]水的可壓縮性時，有的文獻采用聲波單元，有的則采用液體單元，所得結(jié)果雖有差別，但總的趨勢是考慮庫水壓縮性和水庫邊界吸收作用后，大壩的變形和應(yīng)力均有所降低。

圖6 各種模型計算重力壩壩面動水壓力的頻響曲線的對比（L=100m）Fig. 6 Comparison of the calculated frequency response spectra of the hydrodynamic pressure acting on the upstream face of gravity dam using various models

圖7 擴展水庫的重力壩動水壓力頻響曲線（L=1800m）Fig. 7 Frequency response spectra of the hydrodynamic pressure acting on gravity dam for extended reservoir（L=1800m）

我們的看法是在壩—庫水動力相互作用分析中考慮水的可壓縮性與水庫邊界的吸收作用是合理的。在高壩動水壓力分析中，應(yīng)逐步改變采用Westergaard經(jīng)典公式和附加質(zhì)量的計算方法。但希望加強現(xiàn)場觀測，取得必要的實際觀測資料加強檢驗。

3 壩—地基動力相互作用

壩—地基相互作用對大壩的地震響應(yīng)產(chǎn)生重要影響，主要基于以下三方面的因素：

（1）相互作用使壩—庫水—地基系統(tǒng)的振動頻率和振動模態(tài)發(fā)生變化；

（2）系統(tǒng)的振動能量向無限地基進行散發(fā)，產(chǎn)生輻射阻尼；

（3）壩基的地震動輸入較建壩前自由場的地震動發(fā)生了很大的變化。

壩—地基動力相互作用分析仍然是一個比較復(fù)雜的技術(shù)課題。所以在早期的大壩地震響應(yīng)分析中采用了一些簡化的假定，例如采用無質(zhì)量地基的計算模型，略去了輻射阻尼和地震動輸入改變的影響。20世紀80年代以后，壩—地基動力相互作用問題受到了研究者的重視，提出了很多計算模型。其中包括子結(jié)構(gòu)法和基于能量傳遞邊界的直接法等。情況看起來有所改進，但實際上仍然有很多問題有待解決。目前比較普遍的做法是在壩的近場選擇一塊地基與壩一起進行整體分析。例如，在第12屆大壩數(shù)值分析基準研討會上[7]對220m高的拱壩（頂寬約430m，底寬約80m），計算域選為長×寬×高=1000m×1000m×500m。實際的問題是計算域范圍的選擇與計算域截斷邊界采用的能量傳遞模型有關(guān)，目前文獻中的很多傳遞邊界模型屬于局域性的能量傳遞邊界，不考慮節(jié)點間時間與空間相互作用的耦合影響，不少模型近似度比較高，計算可靠性缺乏嚴格的檢驗。還有一個問題是地震動輸入問題，建壩后輸入地震動的變化問題缺乏研究。

目前有關(guān)子結(jié)構(gòu)法或直接法計算大壩與無限地基動力相互作用的計算模型大都建立在均質(zhì)無限地基假設(shè)的基礎(chǔ)上，實際上很多大壩地基的地質(zhì)構(gòu)造都很復(fù)雜，為了適應(yīng)復(fù)雜不均勻地基的動力相互作用分析，我們對圖8所示的幾種典型的復(fù)雜地基，分別提出了相適應(yīng)的計算模型，可供參考。

3.1 層狀不均勻地基

圖8 典型的不均勻地基模型Fig. 8 Typical inhomogeneous foundation model

水平分層地基是一種常見的不均質(zhì)地基模型。我們提出了各向同性和非各向同性層狀不均質(zhì)地基動力剛度的計算方法[8][9]，建立了層狀地基中力和變形的關(guān)系，可適用于結(jié)構(gòu)和地基的動力相互作用分析。對重力壩的分析可直接應(yīng)用。對任意形狀河谷中的拱壩（見圖8a），則可參考文獻[10]建立河谷界面點的格林函數(shù)，求得相應(yīng)的動力剛度矩陣，從而可進行拱壩與地基的動力相互作用分析。

3.2 分塊地基

對于任意形狀河谷的分塊地基（見圖8b），可以方便地采用比例邊界有限元法求解，同時還可以考慮地基剛度（彈性模量E或剪切模量G）和質(zhì)量密度（ρ）隨徑向深度按指數(shù)函數(shù)規(guī)律而變化的情況[11]。應(yīng)用于拱壩的相互作用分析，則需將拱壩地基按無限椎體模型加以一定的改造[12]，如圖9所示。

圖9 拱壩與無限地基動力相互作用的計算模型Fig.9 Computational model for dynamic interaction analysis between arch dam and unbounded foundation

3.3 近場不均勻性不規(guī)則分布的地基

如圖8c所示，壩基近場含不規(guī)則的軟弱區(qū)域，如捕虜體時，可以采用阻尼影響抽取法進行分析。方法的基本思想?yún)⒁妶D10。將壩和地基劃分為兩個子結(jié)構(gòu)，對含不規(guī)則分布的不均勻性地基人工加入高阻尼。高阻尼的存在導(dǎo)致振動波的波幅將很快衰減，從而可以認為含高阻尼有限域的動剛度逼近無限域的動剛度，再進一步設(shè)法將高阻尼影響消去就得到未引入高阻尼時的無限域地基動剛度[11]。

基本計算式表示見式（9）:

式（9）中[S∞（ω）]代表無限域地基的動剛度；[Sζ（ω）]代表引入高阻尼后有限域的動剛度，由于波動能量衰減較快，可以認為其逼近無限域的動剛度；式（9）右端第二項代表對阻尼影響的改正，ζ為引入的人工阻尼。式（9）是頻率域的表達式，為數(shù)值分析方便，要求轉(zhuǎn)換為時域的表達式。為此，對式（9）進行Fourier變換，得出時域位移單位脈沖響應(yīng)函數(shù)的表達式（10）：

圖10 阻尼影響抽取法計算模型Fig. 10 Computational model for damping solvent extraction method

據(jù)此可以得到大壩和地基接觸面上動力相互作用力的表達式。

數(shù)值分析中遇到的問題是應(yīng)該選擇多大的人工阻尼與相適應(yīng)的計算域大小。過大的計算域所需的計算工作量很大，計算有困難，較小的計算域要求引入較高阻尼，阻尼過高相互作用力的特性將發(fā)生變化，對計算結(jié)果的精度和可靠性發(fā)生影響。為克服這一困難，我們提出了阻尼影響逐步抽取的方法[13]，即將阻尼分n次引入和抽取，ζ=nΔζ。每次引入的阻尼量為Δζ，這相當于將式（10）改寫為：

據(jù)此，大壩—地基接觸面上的動力相互作用力{Rb（t）}可表示為 ：

將式（11）代入上式得計算關(guān)系式：

其中：

實際應(yīng)用時，為了減輕計算工作量，{RbL（t）}將不由卷積式（14）推求，而可由含阻尼Δζ的運動方程方便地解出。

數(shù)值算例表明，將高阻尼分多次引入和抽取，既保證了計算精度，又提高了計算效率。

4 大壩抗震分析的精細化數(shù)值模型

大壩的地震應(yīng)力目前仍然是大壩抗震安全評價的重要依據(jù)。有限元方法（FEM）由于其對本構(gòu)方程、邊界條件、幾何形狀和材料特性的靈活適應(yīng)性，在壩工設(shè)計中獲得廣泛應(yīng)用。但FEM用于壩工設(shè)計也有所不足，F(xiàn)EM采用低階插值函數(shù)，例如，我國大壩分析軟件普遍采用線性單元，單元連續(xù)性差，應(yīng)力計算值對網(wǎng)格的依賴性強。按不同離散的單元網(wǎng)格形狀和單元尺寸大小進行分析，將得出不同的單元應(yīng)力計算值。從而應(yīng)力計算精度低、可靠性差。

我國《混凝土拱壩設(shè)計規(guī)范》和《水工建筑物抗震設(shè)計規(guī)范》都規(guī)定在進行有限元法分析時，還需要同時進行拱梁模態(tài)法分析，希望從多方面來彌補這一缺點。國際大壩委員會舉辦的第12屆大壩數(shù)值分析的基準研討會規(guī)定，所有論文有關(guān)的拱壩壩體、水庫、地基單元均采用二階連續(xù)性單元（20節(jié)點四邊形六面體或12節(jié)點三角形五面體），目的在于提高單元的連續(xù)性，以加強成果的可比性。但有限元法采用高階單元，其計算復(fù)雜性和計算工作量將相應(yīng)增加。

為了提高應(yīng)力計算的精度，我們嘗試過多種途徑，等幾何分析（Isogeometric Analysis）是其中之一，這是近期發(fā)展的比較熱門的一種數(shù)值分析方法，可以精確描述形體的幾何形狀，實現(xiàn)自動化網(wǎng)格細分，便于和計算機輔助設(shè)計相結(jié)合，達到較高的應(yīng)力計算精度。比較起來，我們推薦的比例邊界有限元方法SBFEM可以作為提高應(yīng)力計算精確性的基本方法。這是在FEM基礎(chǔ)上發(fā)展起來的一種有廣闊應(yīng)用前景的新型數(shù)值計算方法，可以和FEM實現(xiàn)無縫連接，易于為廣大工程技術(shù)人員所接受。這種方法有很多特點和優(yōu)點。首先SBFEM只需進行邊界離散，使問題降階一維，可節(jié)省很多計算工作量。其次，二維和三維單元的插值函數(shù)可由二維或三維單元邊界上的一維插值函數(shù)所構(gòu)成，所以，構(gòu)造簡單，升階方便，單元間和單元內(nèi)可以實現(xiàn)高階連續(xù)性，這就為單元的應(yīng)力計算達到理想精確度提供了有利的條件。單元邊線上的節(jié)點數(shù)，也就是單元的自由度數(shù)可以自由選擇?？梢愿鶕?jù)需要采用含任意邊數(shù)的多邊形單元。關(guān)鍵是增加單元的節(jié)點數(shù)和提升插值函數(shù)的階次（連續(xù)性），并不增加復(fù)雜性，工作量的增加也很有限，這都是提高單元連續(xù)性，提高應(yīng)力計算精度的有利條件。這也使粗細網(wǎng)格過渡不增加任何復(fù)雜性。

關(guān)于粗細網(wǎng)格過渡問題還可以進一步說明如下。塊體單元與板單元的連接，對于有限元方法來說，是個非常復(fù)雜的問題，因為塊體單元每個節(jié)點只有三個平動自由度，即x、y、z三個方向的變形u、v、w，而板單元的節(jié)點自由度多于3個，除平動外，還要包括節(jié)點的轉(zhuǎn)動自由度，這就增加了板與塊體單元連接的難度和復(fù)雜性。而采用SBFEM，處理就比較方便。應(yīng)用SBFEM進行板結(jié)構(gòu)的分析，每個節(jié)點只需要三個平動自由度u、v、w就可達到很高的計算精度，所以，板與塊體單元的連接，無須增加節(jié)點自由度，這使板與塊體單元的連接非常方便。再者，由于塊體單元可以任意增加邊界面上的節(jié)點（或者自由度），所以，無論板的厚度如何，無論塊體單元本身的尺寸如何，都可直接連接，而不改變板或者塊體單元原來網(wǎng)格的拓撲結(jié)構(gòu)。也就是說，采用SBFEM進行分析，不僅粗、細網(wǎng)格過渡很方便，而且板與塊體的過渡也很方便。關(guān)于殼體的分析，目前，SBFEM還沒有很多經(jīng)驗，這方面的研究正在開展，暫時可將殼近似地看作由板結(jié)構(gòu)組合而成。

SBFEM的主要特點還在于這是半解析半數(shù)值的計算方法，在徑向具有解析解，環(huán)向可達到FEM的離散精度。由于SBFEM單元的高階離散比FEM簡便得多，所以，SBFEM環(huán)向的計算精度也高于FEM。SBFEM的半解析特點，使其對于求解無限域和奇異性（應(yīng)力集中）問題具有特殊的優(yōu)越性。前面所提到的大壩和庫水的動力相互作用問題，大壩和無限地基的動力相互作用問題都屬于無限域問題。對于拱壩和庫水的動力相互作用分析，有限元方法需要離散三維水庫，而SBFEM只需離散壩面。拱壩和無限地基的動力相互作用分析，對于均質(zhì)無限地基，SBFEM可以只在壩和地基的接觸面節(jié)點處進行離散，而FEM則需離散足夠大范圍的地基。

SBFEM在徑向具有解析解，這為求解應(yīng)力集中問題創(chuàng)造了有利條件。例如，斷裂問題裂縫尖端的應(yīng)力強度因子可以解析求出。將這一特點與多邊形單元相結(jié)合，進行大壩斷裂分析時，可以很方便地跟蹤裂縫擴展路徑，而不必像FEM那樣需在裂尖進行單元細分以便計算應(yīng)力強度因子。

SBFEM也有一定的局限性和不足之處，亦即計算域材料特性可以是非各向同性的，但必須是均質(zhì)的，如庫水，均質(zhì)無限地基等。此外，進行非線性分析還有一定困難，由于SBFEM與FEM可以進行無縫連接，可以對易產(chǎn)生非線性響應(yīng)的部分采用FEM分析，而其他單元采用SBFEM分析，這樣可以提高整個問題的計算精度。SBFEM正處于發(fā)展階段，很多問題都可以進行探索。

下面通過幾個實際算例使讀者進一步了解SBFEM。

4.1 粗細網(wǎng)格過渡和不同結(jié)點單元之間的過渡[14]

以薄壁圓筒承受均勻內(nèi)壓作用為例進行分析，問題具有解析解。由于對稱性可以離散圓筒的1/4。圖11表示SBFEM（采用多邊形單元）和FEM兩種方法的離散網(wǎng)格，圖12表示兩種方法與解析解相比的相對誤差與收斂性的效果。由圖可見采用SBFEM節(jié)點自由度數(shù)顯著減少，同時網(wǎng)格細化后，收斂效率高。圖中SBP代表比例邊界有限元多邊形網(wǎng)格，SBP-NECB代表比例邊界有限元與IGA的曲線邊界相結(jié)合。

圖11 厚壁圓筒承受均勻內(nèi)壓作用的網(wǎng)格劃分Fig. 11 Mesh discretization for thick wall cylinder loaded by uniform inner pressure

4.2 靜力問題—壩踵奇異應(yīng)力場分析

采用SBFEM進行奇異應(yīng)力分析時，要求將相似中心置于奇異應(yīng)力點（應(yīng)力集中點）處，因此，對計算方法需要做一定的處理[15]。首先按SBFEM求出應(yīng)力奇異點所在單元周邊的邊界條件，再將相似中心置于應(yīng)力奇異點，求出奇異應(yīng)力場。現(xiàn)選擇Koyna重力壩在自重和水壓力作用下壩踵處的奇異應(yīng)力場分析為例加以說明（見圖13）。

壩的剖面和FEM網(wǎng)格剖分如圖13所示。SBFEM采用多邊形網(wǎng)格，壩踵區(qū)單元為ABC，單元各邊的節(jié)點數(shù)可以調(diào)節(jié)（見圖13，圖14）；有限元網(wǎng)格，在壩踵區(qū)附近進行加密。為了研究壩踵區(qū)的奇異應(yīng)力場，采用了三種類型的網(wǎng)格：①SBFEM網(wǎng)格，節(jié)點分布密度如圖 14中的（a）、（b）、（c）所示；②SBFEM網(wǎng)格，相似中心位于應(yīng)力奇異點A，并建立AC、AB邊的應(yīng)力或位移邊界條件，如圖14中的（d）、（e）所示；③FEM網(wǎng)格，加密情況如圖14中的（f）、（g）、（h）所示。各種網(wǎng)格計算的豎直應(yīng)力σy沿壩基接觸面AB的應(yīng)力分布如圖15所示。由圖可見，經(jīng)過改進處理的SBFEM方法（圖中的Case 2和Case 1）可以比較少的自由度獲得比較準確的奇異應(yīng)力場的計算結(jié)果。此外，還進一步計算了廊道周邊的第一主應(yīng)力分布，如圖16所示。有限元不同的網(wǎng)格密度得出的應(yīng)力分布差別很大，無限加密后才能得出準確的結(jié)果；但采用比例邊界有限元多邊形網(wǎng)格時，廊道凹角處的應(yīng)力分布可以接近于有限元無限加密后的結(jié)果。

圖12 SBFEM與FEM計算精度與網(wǎng)格節(jié)點數(shù)的對比Fig. 12 Comparsion of computational accuracy and number of mesh nodes needed by using FEM and SBFEM

圖13 柯伊那壩的SBFEM和FEM離散網(wǎng)格Fig. 13 Mesh discretization for Koyna dam analysis using SBFEM and FEM

圖14 壩踵區(qū)奇異應(yīng)力場分析采用的各種計算網(wǎng)格Fig. 14 Various mesh discretization in the dam heel region for singular stress field analysis

圖15 各種方法壩基接觸面奇異應(yīng)力場計算結(jié)果比較Fig. 15 Comparison of the calculated results of singular stress field along the dam-foundation interfac obtained by various models

圖16 各種模型廊道周邊應(yīng)力分布計算結(jié)果的比較Fig. 16 Comparison of the calculated results of stress distribution along the surround of the gallery obtained by various models

4.3 動力問題—壩踵奇異應(yīng)力場分析

動力問題的奇異應(yīng)力場一般也可以采用靜力問題類似的方法進行分析。但需要指出，動力問題相對要復(fù)雜得多。采用網(wǎng)格細分等加密方法以后，奇異點附近的應(yīng)力分布呈現(xiàn)波動性收斂。

仍然選擇柯依那壩在柯伊那地震作用下（順河向和豎向地震加速度分別為0.474g和0.316g）的動力響應(yīng)進行分析，考慮大壩—庫水—地基動力相互作用的影響。采用基于等幾何分析IGA（Isogeometric Analysis）的網(wǎng)格，這是目前各種數(shù)值計算方法中相對最為準確的方法，網(wǎng)格可以自動加密。各種網(wǎng)格的計算自由度如表1所示。IGA網(wǎng)格加密至自由度數(shù)達到45030后，壩踵附近奇異應(yīng)力場仍然呈現(xiàn)波動的特點，如圖18所示。采用我們推薦的經(jīng)過處理后的SBFEM多邊形方法，其應(yīng)力收斂趨勢則是一致和協(xié)調(diào)的，同時曲線比較光滑。

表1 各種計算網(wǎng)格的自由度數(shù)（DOFs）和CPU時間比較Tab.1 DOFs and CPU times for various mesh discretizations

圖17 大壩—庫水—地基系統(tǒng)動力響應(yīng)分析的計算網(wǎng)格Fig. 17 Computational mesh for dynamic response analysis of dam-reservoir-foundation system

圖18 壩—基接觸面第一主應(yīng)力峰值分布Fig. 18 Distribution of the peak first principle stress along the dam-foundation interface

4.4 裂縫開合縫內(nèi)水壓力變化對大壩裂縫擴展的影響

大壩上游面裂縫在地震作用下發(fā)生開合，縫內(nèi)水壓力隨之發(fā)生變化，這將對裂縫的擴展發(fā)生影響。采用SBFEM處理縫內(nèi)水壓力非常方便[16]。我們根據(jù)柯依那壩的震害觀測，假設(shè)該壩上游面60m深處出現(xiàn)一條長1m左右的裂縫，研究其進一步的擴展。裂縫內(nèi)水壓隨裂紋嘴開合的變化采用Javamardi等提出的模型[17]進行模擬。地震作用下裂縫面快速地開合導(dǎo)致縫內(nèi)水體不能及時排出，從而形成較大的水壓力。這表明裂縫處在張開、靜止和閉合狀態(tài)下，縫內(nèi)水壓分布顯著不同，張開狀態(tài)下的縫內(nèi)水壓力的強度僅達到閉合狀態(tài)下的8.2%左右。我們的計算結(jié)果與Pekau的積分計算結(jié)果[18]吻合得較好（見圖19）?？乱滥菈卧诘卣鹱饔孟聣巍畮臁鼗到y(tǒng)的計算模型如圖20所示。分析考慮了三種計算模型：①模型1，不考慮大壩—庫水—地基的動力相互作用，采用慣性力輸入，按Westergaard經(jīng)典公式計算附加質(zhì)量，同時不考慮縫內(nèi)水壓力的作用；②模型2，考慮大壩—庫水—地基的動力相互作用，但不考慮縫內(nèi)水壓力的作用；③模型3，考慮大壩—庫水—地基的動力相互作用，同時又考慮縫內(nèi)水壓隨裂縫的張開、閉合的變化。三種模型計算得到的裂縫擴展情況各不相同，參見圖21。第三種模型縫內(nèi)動水壓力隨裂縫的張合發(fā)生劇烈而復(fù)雜的變化，產(chǎn)生的裂縫擴展路徑也最長。這種影響不容忽視。

4.5 多層疊合板的分析

圖19 裂縫嘴張開、閉合狀態(tài)下縫內(nèi)水壓力的分布（我們與Pekau計算結(jié)果的對比）Fig.19 Hydrodynamic pressure distribution inside the crack during the opening and closing state of crack mouth（comparison of the results obtained by the proposed approach and that by Pekau）

圖20 Koyna壩比例邊界有限元多邊形單元計算模型Fig. 20 Discretization of scaled boundary polygon elements for Koyna dam analysis

比例邊界有限元法可以很方便地進行薄板和中厚板的結(jié)構(gòu)分析，每結(jié)點3個平動自由度u、v、w，不必增加結(jié)點的轉(zhuǎn)動自由度。表2所示為我們采用SBFEM所進行的3層疊合板（纖維走向0°/90°/0°）分析結(jié)果[19]與文獻采用高階位移（含3次多項式表示的平動與轉(zhuǎn)角位移）表示的剪切變形理論計算結(jié)果[20]，以及理論解的對比。算例為簡支方形板，承受著按正弦波形式分布的荷載作用。表中分別給出不同板的邊長與板厚比例l/h的情況下規(guī)格化的板中心點位移與應(yīng)力的計算結(jié)果。從表中數(shù)據(jù)可見，采用SBFEM分析可以達到很高的計算精度。由于板的分析中每個節(jié)點可以只采用3個平動自由度，這給大壩結(jié)構(gòu)分析中板與塊體的連接（例如薄閘墩與壩體的連接）提供了很大的方便，不必在連接部位進行網(wǎng)格加密和過渡。

圖21 不同計算模型預(yù)測的Koyna大壩裂縫擴展路徑Fig. 21 Predicted crack path of Koyna dam obtained by various computational models

圖22 多層疊合板分析示意圖Fig. 22 Schematic drawing of laminated composite plates

表2 方形疊合板（0°/90°/0°）的位移和應(yīng)力對比Tab. 2 Nondimensionalized defections and stresses in a three-layer（0°/90°/0°）simple supported square liminate plate subjected to sinusoidal transverse load

5 大壩的超載潛力與抗震安全性評價

大壩在極端地震作用下的非線性動力響應(yīng)與損傷破壞的發(fā)展受到各國研究者的重視，在這方面進行深入研究的興趣也與日俱增?？梢哉J為，這個問題目前還正處在發(fā)展階段。大連理工大學(xué)也進行了這方面的探索，談一點初步看法。

對于大壩的非線性地震響應(yīng)目前可以采用彈塑性分析，斷裂分析等各種數(shù)值計算方法，但準確預(yù)測大壩的破壞模態(tài)以及超載潛力的定量化仍然存在困難。據(jù)了解，美國墾務(wù)局2006年11月曾提出過一份報告《大壩抗震安全急待解決的問題》（草案），其中指出，計算壩的非線性響應(yīng)在很大程度上受所采用的阻尼和材料強度值的影響。雖然數(shù)值計算能力可以很先進，但在非線性階段以及接近破壞階段的輸入變量，目前還受到不確定因素的困擾。有必要進行大比尺的模型試驗來檢驗現(xiàn)行的結(jié)構(gòu)分析模型和方法，包括：初始損傷判定，適宜的阻尼水平，河谷不均勻的地震動輸入，以及幾何非線性（橫縫）和材料非線性引起的壩的非線性響應(yīng)等。

美國墾務(wù)局的報告指出大壩地震破壞模態(tài)與超載潛力預(yù)測存在的問題在于缺乏實際檢驗。一方面經(jīng)歷過高強地震作用的大壩為數(shù)甚少，另一方面獲得的實際觀測資料也很不完整。報告建議通過大比尺的模型試驗來解決問題。認為大比尺的模型可以模擬拱壩的抗剪鍵槽，未膠結(jié)的澆筑面和伸縮縫等。不過我們認為，振動臺上的大比尺模型試驗實際也具有局限性。模型試驗在20世紀對混凝土大壩的建設(shè)曾發(fā)揮過很重要的作用，但那時試驗主要模擬的是大壩的彈性應(yīng)力，要模擬大壩的破壞模態(tài)，有三個方面的問題還有待解決。首先，是仿真模型材料的選擇，不僅要求材料各階段的阻尼相似，特別是接近破壞時的動態(tài)特性和速率影響的模擬。其次，地基和水庫無限性，水庫邊界對波的吸收作用，壩基地震動輸入不均勻的影響是很難達到相似的，更不用說由于水的壓縮性所引起的動水壓力相位差影響的模擬。第三是荷載的相似性。模型水庫中的水一般采用天然水進行模擬。由于模型小水壓力低，為了滿足相似要求，所施加的地震力也不可能很大。要達到模型壩破壞，必須不斷加大所施加的地震力，但是靜水壓力不可能也隨之加大。最終達到破壞時，地震慣性力和靜水壓力的比例與實際情況是不可能相符的。除非采用人工設(shè)計的靜水壓力的施加裝置。

這樣看來，在現(xiàn)有的條件下，對大壩破壞模態(tài)和超載潛力的評估，只能通過彈塑性分析、斷裂分析、其他類型的地震損傷破壞的數(shù)值模擬方法（例如，我們所提出的細觀仿真模擬方法[21]）、仿真材料動力模型試驗、強地震災(zāi)區(qū)大壩震損分析和安全評價等多種途徑的綜合判斷加以解決。

值得重視的發(fā)展動向是美國墾務(wù)局報告中所指出的目前流行的新思路—風險分析。這是從不確定性的觀點出發(fā)，綜合考慮地震發(fā)生的概率，大壩發(fā)生各級損傷的概率，材料強度和材料動態(tài)特性的概率，建立大壩發(fā)生各級損傷的易損性曲線，以便對大壩發(fā)生各級損傷所可能帶來的人員傷亡和經(jīng)濟損失做出評估。這方面我們也曾經(jīng)做過一些工作[22]，評價了大壩發(fā)生各級損傷的概率，但無法評估各級損傷所產(chǎn)生的損失后果，因為缺乏有關(guān)資料。進行這種評估需要在大量調(diào)查研究的基礎(chǔ)上進行。

[1] CHOPRA A.K.， Earthquake Analysis of Arch Dams: Factors to be Considered.[C] // The 14th World Conference on Earthquake Engineering， 2008， Beijing， China.

[2] CHOPRA A.K.， Earthquake Analysis， Design and Safety Evaluation of Concrete Arch Dams[C] // Proceedings of 10th WCEE， Vol.11， 6763-6772.

[3] HALL J.F.， Chopra A.K.， Dynamic Analysis of Arch Dams Including Hydrodynamic Effects， J. ENG. MECH. Div.ASCE 109， 1983，149-167.

[4] TAN H， CHOPRA A.K.， Earthquake Analysis of Arch Dams Including Dam-Water-Foundation Rock Interaction[J]. Earthquake Engineering and Structural Dynamics，1995，24: 1453-1474.

[5] LIN G.， WANG Y.， HU Z. An Efficient Approach for Frequency-Domain and Time-Domain Hydrodynamics Analysis of Dam-Reservoir Systems[J]. Earthquake Engineering and Structural Dynamics，2012，41:1725-1749.

[6] BOUAANANI N， MIGUEL B.， A new Formulation and Error Analysis for Vibrating Dam-reservoir Systems with Upstream Transmitting Boundary Conditions[J]. Journal of Sound and Vibration， 2010， 329 :1924-1953.

[7] ICOLD Proceedings， 12th International Benchmark Workshop on Numerical Analysis of Dams， Oct 2013， Graz-AUSTRIA，Edited by Zenz G and Goldgruber M.

[8] 林皋，韓澤軍，李建波. 層狀地基任意形狀剛性基礎(chǔ)動力響應(yīng)求解[J].力學(xué)學(xué)報，2012，44（6）：1016-1027.LIN G，HAN Zejun.，LI Jianbo， Solution of the Dynamic Response of Rigid Foundation of Arbitrary Shape on Multilayered Soil[J]. Chinese Journal of Theorectical and Applied Mechanics，2012，44（6）：1016-2017.

[9] LIN G.， HAN Z.，LI J.， General Formulation and Solution Procedure for Harmonic Response of Rigid Foundation on Isotropic as well as Anisotropic Multilayered Half-Space[J]. Soil Dynamics and Earthquake Engineering，2015，70：48-59.

[10] 林皋. 地下結(jié)構(gòu)地震響應(yīng)的計算模型[J].力學(xué)學(xué)報（將發(fā)表）.LIN Gao. Computational Model for Seismic Response Analysis of Underground Structures[J]. Chinese Journal of Theoretical and Applied Mechanics.（to be publised）.

[11] WOLF J.P.， SONG C.， Finite-Element Modelling of Unbounded Media， 1996， John Wiley & Sons.

[12] LIN G.， DU J.G.， HU Z.Q.， Earthquake Analysis of Arch and Gravity Dams Including the Effects of Foundation Inhomogeneity[J]. Frontiers of Architecture and Civil Engineering in China， 2007， 1:41-50.

[13] YIN X.， LI J.， WU C.， LIN G.， ANSYS Implementation of Damping Solvent Stepwise Extraction Method for Nonlinear Seismic Analysis of Large 3-D Structures[J]. Soil Dynamics and Earthquake Engineering， 2013， 44:139-152.

[14] 林皋，龐林.大壩結(jié)構(gòu)靜動力分析的精細化模型[J].地震研究，2016，39（1）：1-9.LIN Gao，PANG Lin. Refined Model for Atatic and Dynamic Analysis of Dam Structures[J]. Journal of Seismological Research，2016，39（1）:1-9.

[15] PANG L.， LIN G.， HU Z.， A Refined Global-local Approach for Evaluation of Singular Stress Field Based on Scaled Boundary Finite Element Method， Acta Mechanica Solida Sinica， 2017.

[16] 龐林，林皋.縫水壓力和裂縫面接觸條件對重力壩裂縫擴展的影響[J].計算力學(xué)學(xué)報.（將發(fā)表）PANG Lin，Lin Gao. Crack Development of Gravity Dam Influenced by the Water Pressure Inside the Crack and the Contact State of Cracked Faces[J]. Chinese Journal of Computational Mechanics（to be published）

[17] JAVANMARDI F， LEGER P， TINAWI R， Seismic Water Pressure in Cracked Concrete Gravity Dams: Experimental Study and Theoretical Modeling [J]， Journal of Structural Engineering，2005 ， 131（1）: 139-150.

[18] PEKAU O A， ZHU X. Effect of Seismic Uplift Pressure on the Behavior of Concrete Gravity Dams with a Penetrated Crack[J].Journal of Engineering Mechanics， 2008， 134（11）:991-999.

[19] ZHANG P C， LIN G， Analysis of Laminated Composite Plates Base on the Scaled Boundary Finite Element Method.（to be published）

[20] KANT T， SWAMINATHAN K. Analytical Solutions for Static Analysis of Laminated Composite and Sandwich Plates Based on a Higher-order Refined Theory[J]. Composite Structures， 2002，56: 329-344.

[21] ZHONG H， LIN G， LI X， LI J B. Seismic Failure Modeling of Concrete Dams Considering Heterogeneity of Concrete[J]. Soil Dynamics and Earthquake Engineering， 2011， 31（12）:1678-1689.

[22] Chinese National Committee on Large Dams， Seismic Safety of Dams in China，2014.

Seismic Analysis and Safety Evaluation of Large Dams

LIN Gao
（Dalian University of Technology，Dalian 116024， China）

Based on the research work recently carried out at Dalian University of Technology， some problems which are important for the seismic analysis and safety evaluation of large dams are studied and discussed. The main contents include dynamic damreservoir interaction， dynamic dam-foundation interaction， the refined numerical method for stress analysis and the prediction of the overload capacity and safety evaluation of large dams.

scaled boundary finite element method;earthquake resistant behaviour of large dams;dynamic dam-reservoir interaction; dynamic dam-foundation interaction;high accuracy stress analysis

TV642.1

A

570.25

10.3969/j.issn.2096-093X.2017.02.002

國家重點研發(fā)計劃復(fù)雜工程力學(xué)高性能應(yīng)用軟件系統(tǒng)研制（2016YFB0201001）。

2017-03-19

2017-03-21

林 皋（1929—），大連理工大學(xué)教授，中國科學(xué)院院士，主要研究方向：大壩抗震、核電結(jié)構(gòu)抗震等方面。在中、外學(xué)術(shù)刊物上發(fā)表研究論文600余篇。