基于SBFEM的豎向地震重力壩動水壓力算法研究

王 毅， 林 皋， 胡志強

(大連理工大學(xué) 建設(shè)工程學(xué)部，遼寧 大連 116024)

壩水相互作用問題一直以來都是研究的重點，從20世紀(jì)30年代開始，許多學(xué)者都致力于對壩水動力相互作用的研究。Westergaard[1]最早推導(dǎo)了動水壓力的經(jīng)典附加質(zhì)量理論公式，為地震作用下動水壓力的研究做出了開創(chuàng)性的工作。20世紀(jì)50年代，有限元方法的發(fā)展，對壩水相互作用問題的研究起到了極大的推動作用，其中最為著名的就是以Chopra等[2-5]在對重力壩以及拱壩的壩水動力相互作用進行了比較完善的研究。隨后又有學(xué)者通過邊界元[6-7]或者邊界元與有限元耦合的方法[8]來求解壩水相互作用問題。Bouaanani 等[9-10]采用解析的方法得到了能夠考慮庫水的壓縮性以及庫底吸收的頻域重力壩動水壓力求解公式。

隨著數(shù)值計算方法的不斷發(fā)展，一種新的半解析數(shù)值方法-比例邊界有限元方法(SBFEM)脫穎而出，由于其在處理無限域波動問題中存在的極大的優(yōu)勢，大大促進了壩水相互作用問題的研究。林皋等[11-15]最早開展了SBFEM對壩水相互作用的研究，全面考慮庫水可壓縮性、水庫邊界吸收等因素以及順河向、橫河向與豎向地震動作用對拱壩、重力壩的水動力相互作用的影響。SBFEM使很多復(fù)雜的壩水相互作用問題的計算在很大程度上得以簡化，并從頻域分析發(fā)展到時域分析，均取得了較好的效果。王翔等[16-18]采用比例邊界有限元方法推導(dǎo)了二維高階雙漸進透射邊界條件，并應(yīng)用在求解大壩庫水相互作用中。李上明[19]采用比例邊界有限元方法與有限元相結(jié)合的方法分析了瞬態(tài)壩庫耦合問題。

先前的理論大多忽略了豎向地震動對重力壩動水壓力的影響，往往只是考慮了順河向的地震動激勵對其的影響，Chakrabarti等[20]給出了豎向地震動的影響，但是其僅僅是考慮了庫水的可壓縮性。著名學(xué)者Chopra[21]對Westergaard公式進行了適用性討論并加以改進，而且還推導(dǎo)了豎向地震動激勵情況下壩面動水壓力的解析公式，討論了豎向地震動激勵對壩面動水壓力的影響。

本文在比例邊界有限元的體系下，重點討論豎向地震動激勵對動水壓力響應(yīng)的影響，不僅能夠考慮庫水的可壓縮性還能夠模擬庫底的吸收作用，而且實現(xiàn)了動水壓力的頻域和時域求解。

1 基于比例邊界有限元方法的動水壓力基本理論

在地震動作用下，考慮水體的可壓縮性，假設(shè)水體為小擾動無黏性液體，整個水域應(yīng)滿足Helmholtz方程

(1)

式中：p表示動水壓力，2為Laplace算子，上標(biāo)“··”表示對時間的二階導(dǎo)數(shù)，c為水中壓力波波速，忽略庫水的可壓縮性，c值就為無窮大，方程(1)退化為標(biāo)準(zhǔn)的Laplace方程，在考慮水的可壓縮性的時候，c值取為1 440 m/s。

圖1 壩水系統(tǒng)示意圖

在壩水交界面上滿足應(yīng)力連續(xù)性邊界條件：

(2)

為了考慮庫底的柔性影響，采用Hall等[22]提出的一維吸收邊界條件：

(3)

在無窮遠處滿足輻射條件：

(4)

本文采用SBFEM算法，由于僅存在壩水交界面的離散邊界，無窮遠處的邊界不需要離散，因此能夠自動滿足無窮遠處的輻射條件，在計算中不需要對此邊界進行特殊處理。

庫水表面，忽略微幅重力波的作用，此時庫水表面的動水壓力為零。

在比例邊界有限元的框架下，可以把相似中心O選擇在壩體的下游無窮遠的地方(如圖1所示O)，此時可以認(rèn)為水庫的上下表面是平行的關(guān)系，從而只需離散壩水交界面上的節(jié)點，從而大大的減少了計算量，得到整個求解域上的關(guān)于動水壓力的積分弱形式。

(5)

通過對積分弱形式的進一步的處理，可以得到在比例邊界有限元框架下的控制方程和邊界條件，積分域V表示整個水域，積分域S表示水庫的各個邊界面，w表示積分權(quán)函數(shù)。文獻[14]給出了其詳細的推導(dǎo)過程及其求解過程，壩面動水壓力的頻域結(jié)果表示為

(6)

其中的各個矩陣變量的含義以及求解方法都在都在文獻[14]中有詳細的求解方法。

為了便于下文的應(yīng)用，此處需要對頻域動水壓力進行整理，把順河向與豎直向地震動加速度激勵寫成一個列向量，p可以表示為

(7)

(8)

因此，壩面上所承受的壓力在僅受到豎直向地震動激勵時

(9)

2 時域求解

考慮到壩水系統(tǒng)的線性假設(shè)，可以對頻域的計算結(jié)果進行一次傅里葉逆變換(FFT)就能夠得到時域的結(jié)果。具體操作步驟如下：

(10)

3 數(shù)值算例

(11)

式中：w為單位水重，g為重力加速度，H為壩前水深，y為壩體迎水面豎向坐標(biāo)，c為水中波速，本文取為1 440 m/s。

圖2 重力壩計算模型[14]

3.1 頻域結(jié)果

3.1.1 激勵頻率影響

忽略庫底的吸收作用，采用本文方法可以得到壩踵處豎向地震動激勵的動水壓力頻響曲線

如圖3所示，橫坐標(biāo)表示無量綱激勵頻率，其中ω1為水庫的一階自振頻率(ω1=πc/2H)；縱坐標(biāo)表示頻域動水壓力與靜水壓力的無量綱比值的絕對值。結(jié)果表明，采用本文方法與解析解吻合良好，即使是在水庫的前三階自振頻率點附近，仍然能夠得到具有較高精度的結(jié)果。另外，計算結(jié)果表明，對壩踵的動水壓力影響較大的頻率區(qū)間在水庫的一階自振頻率之內(nèi)以及水庫的各階自振頻率附近的相對較小帶寬的頻率區(qū)域。不考慮庫底的吸收作用，在激勵頻率達到水庫的各階自振頻率時，壩面上的動水壓力也將出現(xiàn)峰值，表明壩水系統(tǒng)的無阻尼共振特性。當(dāng)激勵頻率超過水庫的一階自振頻率之后，動水壓力響應(yīng)也急劇下降，尤其是激勵頻率在水庫的兩階自振頻率之間的時，壩水響應(yīng)趨于零，其與解析解表達相吻合。

圖3 壩踵處動水壓力頻響曲線

為了更加直觀的表示豎向動水壓力對整個壩面的作用，圖4給出了壩面整體動水壓力的頻響曲線。橫坐標(biāo)仍為無量綱頻率，縱坐標(biāo)為壩體表面在受到豎向地震動激勵情況下的總體動水壓力，解析解仍為文獻[21]所給出的

(12)

式中:H(Ω)為整個壩面上的動水壓力，求解時需要在迎水面進行積分；F0=ρgH2/2，為壩面上的靜水壓力；Ω=ω/ω1，無量綱激勵頻率。

圖4能夠清晰的表示出整個壩面在豎向動水壓力的作用下的響應(yīng)情況：① 在水庫的各階自振頻率點，動水壓力響應(yīng)存在相同的共振現(xiàn)象；② 動水壓力響應(yīng)的敏感頻率區(qū)間在庫水的二階自振頻率以內(nèi)，與壩踵處的頻響函數(shù)相比，在水庫前兩階自振頻率之間，動水壓力響應(yīng)有了顯著的增大，原因在于壩體迎水面上的其他部位對其的貢獻得到了累加；③ 動水壓力響應(yīng)在水庫二階自振頻率之后急劇減小，逼近于零點，表明整個壩面的響應(yīng)都與壩踵處相差不大，僅在庫水的高階自振頻率點附近一個較小的窄帶寬內(nèi)存在一定的共振響應(yīng)值，同時由于共振峰值的存在，本文計算結(jié)果與解析解之間在此點附近也存在一定的差異。

3.1.2 庫底阻尼對壩水相互作用的影響

圖3和圖4的計算結(jié)果都是在忽略阻尼的情況下得到的理想意義的解。事實上，壩水系統(tǒng)也是存在阻尼的，尤其是庫底的淤沙，對地震波有明顯的吸收作用。本文采用庫底反射系數(shù)的方法來模擬庫底的阻尼，通過選擇不同的反射系數(shù)來表示庫底淤沙的阻尼作用，如圖5所示。隨著反射系數(shù)的降低，庫底的吸收作用愈加明顯，尤其是在水庫的自振頻率點處，峰值被明顯削弱。當(dāng)庫底為完全柔性的情況下，壩水相互作用的頻域特征幾乎完全被改變，在整個頻率范圍之內(nèi)，響應(yīng)值都處于一個較低的水平，尤其是在ω<2ω1的頻率區(qū)間內(nèi)，響應(yīng)幾乎呈線性下降趨勢，完全改變了剛性庫底在水庫的第一階自振頻率處出現(xiàn)峰值的特性，此外，完全柔性庫底仍然保留了壩水相互作用在兩階自振頻率中點處的駐值現(xiàn)象。動水壓力峰值也隨著反射系數(shù)的減低而逐漸降低，而且隨著激勵頻率的增大，壩踵處的動水壓力峰值也存在下降的趨勢，但衰減趨勢并不明顯，這與垂直激勵下關(guān)于動水壓力的加速度脈沖響應(yīng)函數(shù)并不衰減[21]密切相關(guān)。綜上所述，庫底阻尼對壩水系統(tǒng)的動水壓力響應(yīng)起著至關(guān)重要的作用，在實際工程應(yīng)用中可根據(jù)庫底的實際情況對庫區(qū)的阻尼進行估計。

3.1.3 水的可壓縮性的影響

以上算例都是在考慮了庫水的可壓縮性得到的結(jié)果，事實上，在相當(dāng)長的一段時間內(nèi)，都是忽略水的可壓縮性來求解壩水系統(tǒng)的動力響應(yīng)的。不考慮水的可壓縮性，整個水庫的水體控制方程將退化為Laplace方程，從而導(dǎo)致系統(tǒng)的反應(yīng)將與頻率無關(guān)，在頻響曲線上將會是一條平行于頻率軸的直線，其初值是以上各個算例得到的初值，即在ω=0處的響應(yīng)值(如圖5所示)。當(dāng)庫底的反射系數(shù)大于0.25時，水的可壓縮性在庫水一階自振頻率以內(nèi)顯著增大了壩面上的動水壓力，相反在高頻階段，其僅在庫水的自振頻率附近的較小頻率域內(nèi)能夠提高壩面上的動水壓力，而且隨著頻率的不斷增大，這種作用也在逐漸降低。如若庫底的反射系數(shù)小于0.25時，忽略庫水的可壓縮性就在整個頻率范圍內(nèi)較大的估計了在豎向地震動激勵下壩面上的動水壓力。在同時考慮庫水的可壓縮性以及庫底的吸收作用所得到的效果與現(xiàn)行工程中所采用的不可壓縮庫水所得效果不盡相當(dāng)，在庫底條件相對較柔的情況下，即反射系數(shù)小于0.5時，前者所得的計算結(jié)果要明顯小于后者；而當(dāng)庫底相對剛性情況下，即反射系數(shù)大于0.5，庫底淤沙層較薄時，前者計算結(jié)果在水庫的一階自振頻率之內(nèi)要遠遠高于后者。

3.2 時域結(jié)果

以上給出的計算結(jié)果均為頻域結(jié)果，在大壩遭受地震作用下，各個時刻的壩面上的動水壓力往往是我們關(guān)注的重點問題，這就要求能夠在時域中求解壩面上的動水壓力，從而確定在整個地震過程中的最不利時刻。

為了更好地與解析解進行比較，本文選取比較具有代表性的三角形脈沖加速度時稱曲線。脈沖激勵的時長為2 s，峰值為1 m/s2。為了便于比較，本文還給出了在大壩受到水平方向上相同的地震動激勵下，壩面動水壓力的變化情況。

圖6 三角形脈沖加速度

計算結(jié)果表明，水平向結(jié)果圖7和豎向結(jié)果圖8與解析解吻合良好。比較水平向與豎向時域結(jié)果的峰值，假設(shè)大壩在水平與豎向收到相同的地震加速度的激勵，發(fā)現(xiàn)大壩壩踵處的動水壓力接近靜水壓力的8%，而豎向情況下則超過11%。在豎向激勵情況得到的峰值比水平向激勵情況下更大。除了近場地震動之外，一般豎向地震動加速度小于水平相的地震加速度。假設(shè)av=(2/3)ah時，兩者在壩踵處的動水壓力也比較接近。因此，豎向激勵是一個不可忽略的影響，不考慮豎向地震動的影響會明顯低估壩面在地震過程中受到的動水壓力。此外計算結(jié)果表明，動水壓力的時程形式與激勵的形式基本吻合，尤其是在壩體受到順河向激勵的情況下。在壩體受到豎向激勵時候，在整個時程上，動水壓力會出現(xiàn)較小幅度的振蕩，如圖8所示。

4 結(jié) 論

通過對豎向地震動激勵對壩面動水壓力的研究，不難做出以下論點：

(1)豎向地震動激勵對大壩的考慮水庫邊界吸收影響時，豎向動水壓力與水平動水壓力存在一定的相位差，地震響應(yīng)的影響是不可忽略的。頻域結(jié)果和時域結(jié)果均表明豎向地震動引起的壩面動水壓力在本文假設(shè)模型下不亞于同等幅值地震動激勵下的順河向結(jié)果，但由于垂直激勵下動水壓力響應(yīng)與邊界條件、淤沙本構(gòu)等有很大關(guān)系，使得此項響應(yīng)值在實測、實驗和數(shù)值計算結(jié)果上存在一定的差異。因此，在進行大壩-庫水系統(tǒng)的動力分析中應(yīng)該對豎向地震動激勵以足夠的重視。

(2)水的可壓縮性對壩體-庫水系統(tǒng)的影響隨著激勵頻率和庫底反射系數(shù)的變化呈現(xiàn)出不同的特性，采用簡單的不可壓縮庫水假設(shè)在處理壩水相互作用問題上具有明顯的局限性。

(3)庫底的柔性作用對地震波的吸收效果顯著，合理地確定庫底的反射系數(shù)能夠更加真實地反映壩水相互作用問題中淤沙層的吸能效應(yīng)。此外，反射系數(shù)法能夠更加方便地體現(xiàn)庫底的不同的淤沙條件。

(4)采用比例邊界有限元方法對豎向地震動激勵下壩面動水壓力的計算僅需要離散壩面單元，計算量較之有限元方法有大幅度的降低，而且對于無窮遠處的輻射條件能夠較好的模擬，計算結(jié)果精度較高，適合推廣應(yīng)用于三維拱壩的計算中。

參 考 文 獻

[1]Westergaard H M. Water pressures on dams during earthquakes [J]. Transactions (ASCE),1933, 98:418-433.

[2]Chakrabarti P, Chopra A K. Hydrodynamic effects in earthquake response of gravity dams [J]. Journal of the Structural Division, 1974, 100(s16):1211-1224.

[3]Porter C S, Chopra A K. Dynamic analysis of simple arch dams including hydrodynamic interaction [J]. Earthquake Engineering and Structural Dynamics, 1981,9(6):573-597.

[4]Porter C S, Chopra A K. Hydrodynamic effects in dynamic-response of simple arch dams [J]. Earthquake Engineering and Structural Dynamics, 1982,10(3):417-431.

[5]Hall J F, Chopra A K. Dynamic analysis of arch dams including hydrodynamic effects [J]. Journal of Engineering Mechanics (ASCE), 1983,109(1):149-167.

[6]Hanna Y G, Humar J L. Boundary element analysis of fluid domain [J]. Journal of Engineering Mechanics Division (ASCE),1982; 108(2):436-450.

[7]Humar J L, Jablonski A M. Boundary element reservoir model for seismic analysis of gravity dams [J]. Earthquake Engineering and Structural Dynamics, 1988, 16(8):1129-1156.

[8]Estorff O V, Antes H. On FEM-BEM coupling for fluid-structure interaction analyses in the time domain [J]. International Journal for Numerical Method in Engineering, 1991, 31(6):1151-1168.

[9]Bouaanani N, Paultre P, Proulx J. A closed-form formulation for earthquake-induced hydrodynamic pressure on gravity dams [J]. Journal of Sound and Vibration,2003, 261:573-582.

[10]Bouaanani N. Miquel 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.

[11]林 皋, 杜建國. 基于SBFEM的壩-庫水相互作用分析[J]. 大連理工大學(xué)學(xué)報, 2005,45(5):723-729.

LIN Gao, DU Jian-guo. Analysis of dam-reservoir interaction based on SBFEM[J]. Journal of Dalian University of Technology, 2005,45(5):723-729.

[12]Lin G, Wang Y, Hu Z. Hydrodynamic pressure on arch dam and gravity dam including absorption effect of reservoir sediments [J]. IOP Conf. Series：Materials Science and Engineering,2010, 10(1)：012234.

[13]Lin G, Du J, Hu Z. Dynamic dam-reservoir interaction analysis including effect of reservoir boundary absorption [J]. Science in China Series E Technological Sciences, 2007, 50(1);1-10.

[14]Lin G, Wang Y, Hu Z Q. An efficient approach for frequency-domain and time-domain hydrodynamic analysis of dam-reservoir systems [J]. Earthquake Engineering and Structure Dynamics, 2012,41(7):1725-1749.

[15]杜建國，林 皋，謝清糧.一種新的求解壩面動水壓力的半解析方法[J]. 振動與沖擊,2008，27(3):31-34,45.

DU Jian-guo, LIN Gao, Xie Qing-liang. A new semi-analytical method for solving hydrodynamic pressure on dams[J]. Journal of Vibration and Shock,2008，27(3):31-34,45.

[16]王 翔, 金 峰. 動水壓力波高階雙漸進時域平面透射邊界I：理論推導(dǎo)[J]. 水利學(xué)報,2011, 42(7)：839-847.

WANG Xiang, JIN Feng. High-order doubly asymptotic time-domain plane transmitting boundary for hydrodynamic pressure I. Theoretical derivation[J]. Journal of Hydraulic Engineering, 2011, 42(7)：839-847.

[17]王 翔, 金 峰. 動水壓力波高階雙漸進時域平面透射邊界II：計算性能[J]. 水利學(xué)報,2011, 42(8)：986-993.

WANG Xiang, JIN Feng. High-order doubly asymptotic time-domain plane transmitting boundary for hydrodynamic pressure II. Computational performance[J]. Journal of Hydraulic Engineering, 2011,42(8):986-993.

[18]Wang X, Jin F, Prempramote S, et al. Time-domain analysis of gravity dam-reservoir interaction using high-order doubly asymptotic open boundary [J]. Computers and Structures, 2011,89(7-8):668-683.

[19]李上明. 基于比例邊界有限元法的壩庫瞬態(tài)耦合分析[J]. 華中科技大學(xué)學(xué)報(自然科學(xué)版),2011,39(9):108-111.

LI Shang-ming. Dam-reservoir transient coupling analysis based on SBFEM[J]. Journal of Huazhong University of Science and Technology(Natural Science Edition),2011,39(9):108-111.

[20]Chakrabarti P, Chopra A K. Hydrodynamic pressures and response of gravity dams to vertical earthquake component [J]. Earthquake Engineering and Structural Dynamics, 1972,1(4):325-335.

[21]Chopra A K. Hydrodynamic pressures on dams during earthquakes [J]. Journal of the Engineering Mechanics Division, Proceeding of the American Society of Civil Engineers, 1967,93(EM6):205-223.

[22]Hall J F, Chopra A K. Two-dimensional dynamic analysis of concrete gravity dams and embankment dams including hydrodynamic effects[J]. Earthquake Engineering and Structural Dynamics, 1982, 10(2), 305-332.