高面板壩地震動(dòng)三維非一致輸入的波函數(shù)組合法①

姚 虞， 劉天云， 張建民

(清華大學(xué)水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084)

高面板壩地震動(dòng)三維非一致輸入的波函數(shù)組合法①

姚 虞， 劉天云， 張建民

(清華大學(xué)水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084)

提出三維問題的波函數(shù)組合法求解不規(guī)則河谷的散射問題，應(yīng)用到高面板壩動(dòng)力有限元計(jì)算的地震動(dòng)輸入中，可以滿足非一致輸入的要求，并回避人工邊界對(duì)外行波的反射問題。重點(diǎn)介紹三維問題的波函數(shù)組合法的理論推導(dǎo)和程序的算例驗(yàn)證，并給出將該方法應(yīng)用到高面板壩地震動(dòng)輸入的基本步驟。

高面板壩； 三維； 非一致； 位移輸入

0 引言

高面板壩的抗震設(shè)計(jì)常常要進(jìn)行動(dòng)力有限元計(jì)算，在地震動(dòng)輸入方面會(huì)面臨兩個(gè)問題：一是如何反映地震波的行波效應(yīng)[1]，另一個(gè)是如何避免人工邊界對(duì)外行波的反射[2]。在目前常用的堆石壩動(dòng)力分析程序中這兩個(gè)問題往往被忽略，代之以均勻輸入和固定邊界。對(duì)于300 m級(jí)的高面板壩，地震動(dòng)作用的非一致性更加明顯，采用簡(jiǎn)單的均勻輸入與實(shí)際有很大差別。

實(shí)現(xiàn)非一致輸入需要求解計(jì)算邊界即不規(guī)則河谷各點(diǎn)的振動(dòng)過程。由于基巖的波阻抗(介質(zhì)密度與波速的乘積)遠(yuǎn)大于壩體，可以通過求解不規(guī)則空心河谷的散射問題得到邊界各點(diǎn)的振動(dòng)過程。

圖1表示瑞利波或任意角度的壓縮波或剪切波入射到半空間不規(guī)則河谷的示意圖。要求解的問題是在任意入射波的作用下河谷邊界任意一點(diǎn)的震動(dòng)過程。若能求解不規(guī)則河谷任意一點(diǎn)的震動(dòng)過程，則該解可作為河谷內(nèi)土石壩地震動(dòng)響應(yīng)有限元求解的輸入。這種輸入方式既考慮了地震動(dòng)作用的非一致性，又回避了人工邊界對(duì)外行波的反射問題。

圖1 半空間不規(guī)則河谷入射波示意圖Fig.1 Incident wave at an irregular valley in a half-space

彈性波的散射問題已有很長(zhǎng)的研究歷史。Sezawa在1927年完成了對(duì)球、圓柱和橢圓柱體在P波入射情況下的散射問題的研究，利用特殊波函數(shù)構(gòu)造了問題的解。之后對(duì)散射問題的研究陸續(xù)展開，并形成了一些方法。Nishimura[3]解決了各向同性介質(zhì)中球體空洞的動(dòng)應(yīng)力集中問題。同一時(shí)期，Pao和Mow等[4]還探討了圓柱體內(nèi)含物的動(dòng)應(yīng)力集中問題。Baron等[5]率先對(duì)彈性介質(zhì)中圓柱形空腔引起的脈沖散射問題進(jìn)行了分析，使用積分變換和波函數(shù)展開法給出了壓縮波散射的解析解；對(duì)于球形內(nèi)含物和空腔的瞬態(tài)效應(yīng)，Mow和Miklowitz進(jìn)行了求解[6-8]。Datta等[9-11]提出采用匹配漸進(jìn)展開的方法分析了半空間內(nèi)含柱形空洞對(duì)P波、SV波和SH波的散射問題，給出了遠(yuǎn)場(chǎng)的漸近表達(dá)式及相應(yīng)的數(shù)值解。Bostrom和Kristensson等[12]研究了半空間內(nèi)球形空洞對(duì)Rayleigh波的散射問題。上世紀(jì)80年代，Pao等[13]運(yùn)用廣義射線法討論了柱體的瞬態(tài)波散射。黎在良和劉殿魁[14]還建立了一個(gè)射線理論用來(lái)解各向異性介質(zhì)中圓柱體對(duì)SH波的散射問題。

近20多年來(lái)，由于計(jì)算機(jī)技術(shù)的廣泛應(yīng)用，用數(shù)值化方法求解彈性動(dòng)力學(xué)問題取得了進(jìn)展。實(shí)踐證明這些方法已成為解決波動(dòng)問題的重要工具。如廖振鵬[15-16]的結(jié)構(gòu)動(dòng)力反應(yīng)的數(shù)值模擬以及近場(chǎng)波動(dòng)的數(shù)值模擬等。另外，對(duì)不規(guī)則的地形，Sanchez-Sesma等[17-19]先后采用邊界元法和波函數(shù)展開法對(duì)二維、三維問題進(jìn)行了分析，Ray 等[20]、劉國(guó)利等[21]分別采用復(fù)變函數(shù)方法、積分變換和波函數(shù)展開相結(jié)合的方法給予分析。

1 方法原理介紹

由彈性力學(xué)中線彈性體位移表示的位移運(yùn)動(dòng)方程(Navier方程，不計(jì)體力)：

(1)

其中：λ、μ為拉梅常數(shù)；u為位移矢量。再考慮Helmholtz定理[22]：

(2)

其中：φ為速度場(chǎng)的標(biāo)量勢(shì)；ψ為矢量勢(shì)。即任意矢量場(chǎng)可分解為一個(gè)無(wú)散場(chǎng)和一個(gè)無(wú)旋場(chǎng)之和。

令 u(1)=φ，u(2)=×ψ，代入Navier方程，得

(3)

故只需滿足

(4)

在頻域內(nèi)，對(duì)于單一頻率的波，式(4)的解可寫為分離變量的形式，即

(5)

將式(5)帶入式(4)，可得

(6)

式(6)稱為亥姆霍茲方程，是頻域內(nèi)的波動(dòng)方程。在給定邊界條件下求解式(6)得到位移勢(shì)函數(shù)之后便可求解位移場(chǎng)。

對(duì)于不規(guī)則河谷散射問題，要求滿足河谷邊界和半空間自由表面面力為零的邊界條件。為得到滿足邊界條件的解，本文提出一種滿足全局邊界條件的波函數(shù)組合法，可以考慮包括瑞利波在內(nèi)的不同地震波的不同角度入射情況。該方法通過組合體波函數(shù)和Lamb問題的解以滿足半空間全局的邊界條件，來(lái)求解不規(guī)則河谷的散射問題，是一種半解析的方法。目前已用該方法編寫了二維和三維河谷的計(jì)算程序，并分別利用半圓形和半球形河谷與理論解作對(duì)比，驗(yàn)證了程序的正確性。

波函數(shù)組合法的基本思路是：

(1) 確定入射波自由場(chǎng)應(yīng)力、位移；

(2) 位移、應(yīng)力波函數(shù)展開；

(3) 確定自由表面波函數(shù)產(chǎn)生的面力；

(4) 在半無(wú)限空間自由表面反向疊加波函數(shù)產(chǎn)生的面力，求解Lamb問題計(jì)算散射邊界次生面力；

(5) 匹配散射邊界處自由場(chǎng)面力、波函數(shù)產(chǎn)生的面力與Lamb解產(chǎn)生的面力，形成散射邊界無(wú)面力狀態(tài)；

(6) 疊加散射邊界處自由場(chǎng)、波函數(shù)和Lamb解產(chǎn)生的位移，得到散射邊界的總位移。

對(duì)于二維波函數(shù)組合法的理論推導(dǎo)及程序驗(yàn)證，作者已在文獻(xiàn)[23]中有所介紹。本文主要介紹三維波函數(shù)組合法。

2 三維波函數(shù)組合法理論推導(dǎo)

球坐標(biāo)下，式(6)有分離變量形式的通解：

(7)

用波函數(shù)表達(dá)三維半空間的應(yīng)力場(chǎng)和位移場(chǎng)的方法與二維的方法基本一致，其表達(dá)式可參見文獻(xiàn)[24]。三維問題的難點(diǎn)在于三維非軸對(duì)稱Lamb問題的求解。

由于波函數(shù)展開法不能滿足水平自由表面的邊界條件，存在殘余正應(yīng)力σ與參與剪應(yīng)力τ。為滿足自由表面條件，提出在表面處施加與殘余面力相反方向的虛擬面力-σ與-τ(圖2)。

圖2 自由表面殘余面力的消除Fig.2 Elimination of residual force on a free surface

波數(shù)域內(nèi)，三維蘭姆問題的控制方程為：

(8)

其中φ和Ψ滿足

(9)

滿足以上兩式的三個(gè)勢(shì)函數(shù)：

(10)

(11)

將三個(gè)邊界應(yīng)力進(jìn)行組合，有

(12)

(13)

解得

(14)

3 算例驗(yàn)證

三維程序利用半球形河谷的散射問題進(jìn)行驗(yàn)證。驗(yàn)證結(jié)果參考文獻(xiàn)[18]，該文獻(xiàn)使用波函數(shù)展開法，將三倍半徑范圍的半空間自由表面納入計(jì)算域邊界中，忽略三倍半徑以外自由表面的邊界條件。文獻(xiàn)中提供了不同頻率P波垂直向上入射的結(jié)果，筆者用本方法計(jì)算出相應(yīng)的結(jié)果以進(jìn)行對(duì)比(圖3)。

圖3 P波垂直入射半球形河谷邊界位移放大系數(shù)Fig.3 Displacement amplification factor of hemispherical valley boundary under vertical incidence of P wave

從圖3中可以看出，采用本文方法計(jì)算的結(jié)果與文獻(xiàn)[18]計(jì)算的結(jié)果接近，從而驗(yàn)證了三維程序的正確性。

4 應(yīng)用方法

地震從遠(yuǎn)處傳到壩址建基面，輸入壩體內(nèi)部，引起壩體的變形和破壞。為模擬地震輸入情況，需要以下條件：(1)入射方位與入射角；(2)入射波形式；(3)地震時(shí)程。

首先，將地震時(shí)程S(t)(一般是間隔0.02 s的離散記錄)通過FFT獲得功率譜曲線W(ω)；

第二，將功率譜曲線W(ω)等間隔離散，如：Δω=0.5，獲得一系列諧波Wi；

第四，對(duì)所有離散頻點(diǎn)進(jìn)行第三步的計(jì)算，可獲得壩址建基面散射場(chǎng)位移譜Ws(x，ω)；

第五，應(yīng)用逆FFT，獲得壩址建基面散射場(chǎng)的位移時(shí)程ws(x，t)。

(15)

式中下標(biāo)s為結(jié)構(gòu)結(jié)點(diǎn)量；b為壩基結(jié)點(diǎn)量；us為結(jié)構(gòu)結(jié)點(diǎn)位移向量。

5 結(jié)論

本文提出的三維波函數(shù)組合法可解決三維半空間內(nèi)任意空心體的散射問題。其基本思路是將地震波轉(zhuǎn)換到頻域之后離散，對(duì)每種頻率的波通過組合波函數(shù)和Lamb問題的解，使得邊界零面力條件得到滿足，再由波函數(shù)和Lamb問題的解得到邊界的頻域內(nèi)的位移，最后轉(zhuǎn)換回時(shí)域，得到邊界的位移時(shí)程。由于堆石體的波阻抗相對(duì)于基巖很小，可將該邊界位移時(shí)程作為壩體動(dòng)力有限元計(jì)算的輸入。

通過對(duì)半球形河谷在P波垂直入射條件下的散射問題進(jìn)行編程計(jì)算，初步驗(yàn)證了理論的正確性。將該方法應(yīng)用到高面板壩地震動(dòng)輸入中，可以考慮地震動(dòng)作用的非一致性，并回避人工邊界對(duì)外行波的散射問題。

References)

[1] 王艷.非一致地震動(dòng)場(chǎng)數(shù)值方法研究及在結(jié)構(gòu)動(dòng)力分析中的應(yīng)用[D].北京：清華大學(xué)，2007.WANG Yan.Research on the Numerical Method for Asynchronous Seismic Wave Motions and its Application in Dynamic Analysis of Structures[D].Beijing：Tsinghua University，2007.(in Chinese)

[2] 李彬.地鐵地下結(jié)構(gòu)抗震理論分析與應(yīng)用研究[D].北京：清華大學(xué)，2005.LI Bin.Theoretical Analysis of Seismic Response of Underground Subway Structures and Its Application[D].Beijing：Tsinghua University，2005.(in Chinese)

[3] Nishimura G.Dynamic Problem of stress Concentration in the Vicinity of a Spherical Matter Included in an Elastic Solid Under Dynamical Force[J].Journal of the Faculty of Engineering，1955，24：101.

[4] Pao Y H，Mow C C.Dynamic Stress Concentration in an Elastic Plate with Rigid Circular Inclusion.No.SR-41[R].The MITRE Corporation，Bedford，Massachusetts，1961.

[5] Baron M L，Mattews A T.Diffaction of a Pressure Wave by a Cylindrical Cavity in an Elastic Medium[J].Journal of Applied Mechanics，1961，28(3)：347-354.

[6] Mow C C.On the Transient Motion of a Rigid Spherical Inclusion in an Elastic Medium and Its Inverse Problem[J].Journal of Applied Mechanics，1966，33(4)：807-813.

[7] Mow C C.Transient Response of a Rigid Spherical Inclusion in an Elastic Medium[J].Journal of Applied Mechanics，1965，32(3)：637-642.

[8] Miklowitz J.Scattering of a Plane Elastic Compressional Pulse by a Cylindrical Cavity[C]// Proceedings of the 11th International Congress of Applied Mechanics.Munich，1966：469-483.

[9] Datta S K，Akily N E.Diffraction of Elastic Waves by Cylindrical Cavity in a Half-space[J].The Journal of the Acoustical Society of America，64(6)：1978：1692-1699.

[10] Datta S K，Shah A H.Truell R，Scattering of SH Waves by Embedded Cavities[J].Wave Motion，1982，4(3)：265-283.

[11] Datta S K.Scattering of Elastic Waves[M]// Nematnasser S, edit.Mechanics Today.New York:Pergamon Press, 1978, 4:149-208.

[12] Bostrom A，Kristensson G.Scattering of a Pulsed Rayleigh Wave by a Spherical Cavity in an Elastic Half Space[J].Wave Motion，1983，5(2)：137-143.

[13] Pao Y H，Ku G C.Application of the Theory of Generalized Rays to Diffraction of Transient Waves by a Cylinder[J].Wave Motion，1983，5(4)：385-398.

[14] 黎在良，劉殿魁.各向異性介質(zhì)中圓柱體對(duì)SH波散射的射線理論[J].地震工程與工程振動(dòng)，1987，7(1)：1-8.LI Zai-liang，LIU Diam-kui.Ray Method of SH Wave Propagation in an Anisotropic Medium with a Circular Cavity[J].Earthquake Enginelring and Engineering Vibration，1987，7(1)：1-8.(in Chinese)

[15] 廖振鵬.近場(chǎng)波動(dòng)的數(shù)值模擬[J].力學(xué)進(jìn)展，1997，26(2)：193-216.LIAO Zhen-peng.Numerical Simulation of Near-field Wave Motion[J].Advances in Mechanics，1997，26(2)：193-216.(in Chinese)

[16] 廖振鵬.結(jié)構(gòu)動(dòng)力反應(yīng)的數(shù)值模擬[M]//中國(guó)科學(xué)技術(shù)前沿(第四卷).北京：高等教育出版社， 2001：427-458.LIAO Zhen-peng.Numerical Simulation of Structure Dynamic Response[M]//Science and Technology at the Frontier in China (Vol.4).Beijing：Higher Education Press，2001：427-458.(in Chinese)

[17] Sanchez-Sesma F J，Herrera I，Aviles J.A Boundary Method for Elastic Wave Diffraction：Application to Scattering of SH Waves by Surface Irregularities[J].Bulletin of The Seismological Society of America，1982，72(2)：473-490.

[18] Sanchez-Sesma F J.Diffraction of Elastic Waves by Three-dimensional Surface Irregularities[J].Bulletin of the Seismological Society of America，1983，73(6)：1621-1636.

[19] Sanchez-Sesma F J，Miguel A B，Ismael H.Surface Motion of Topographical Irregularities for Incident P，SV，and Rayleigh Waves[J].Bulletin of the Seismological Society of America，1985，75(1)：263-269.

[20] Han R PS，Yeh K Y，Liu G L，et al.Scattering of Plane SH-waves by a Cylindrical Canyon of Arbitrary Shape in Anisotropic Media[J].International Journal of Engineering Science,1992,30(12): 1773-1787.

[21] 劉國(guó)利，劉殿魁.位移階躍SH波對(duì)半圓形凹陷地形的散射[J].力學(xué)學(xué)報(bào)，1994：26(1)：70-80.LIU Guo-li，LIU Dian-kui.The Scattering of Displacement Step SH-waves by a Semi-cylindrical Canyon[J].Acta Mechanics Sinica，1994，26(1)：70-80.(in Chinese)

[22] [美]阿肯巴赫，編.徐植信，洪錦如，譯.彈性固體中波的傳播[M].上海：同濟(jì)大學(xué)出版社，1992.Achenbach J D，edit.Translated by XU Zhi-xin，HONG Jin-ru.Wave Propagation in Elastic Solids[M].Shanghai：Tongji Univerity Press，1992.(in Chinese)

[23] 姚虞，劉天云，張建民.求解二維彈性半空間散射問題的波函數(shù)組合法[J].地震工程學(xué)報(bào)，2013，35(1)：109-113.YAO Yu，LIU Tian-yun，ZHANG Jian-min.The Method by the Combination of Eave Functions to Solve the Problem of Scattering in 2D Elastic Half Space[J].China Earthquake Engineering Journal，2013，35(1)：109-113.(in Chinese)

[24] 鮑亦興，毛昭宙.彈性波的衍射與動(dòng)應(yīng)力集中[M].北京：科學(xué)出版社，1993.PAO Yih-hsing，MOW Chao-chow.The Diffraction of Elastic Waves and Dynamic Stress Concentrations[M].Beijing：Science Press，1993.(in Chinese)

Wave Function Combination Method Used for Three-dimensional Non-uniform Input of Earthquake Motion in High Concrete-faced Rockfill Dam

YAO Yu， LIU Tian-yun， ZHANG Jian-min

(StateKeyLaboratoryofHydroscienceandEngineering，TsinghuaUniversity，Beijing100084，China)

The scattering problem of elastic waves is of great engineering significance，especially in earthquake engineering.Frequently，the response of the scatterer is needed to evaluate the safety of structures during earthquakes.The ground motion can be used as dynamic input for the finite-element method calculation，which is a way to realize the non-uniform input.This article puts forward the wave function combination method for a three-dimensional (3D)problem to solve the irregular valley scattering problem.The application of the earthquake input to the dynamic finite-element calculation for a high-face rockfill dam can meet the requirement for non-uniform input and，at the same time，avoid the problem of the reflection of the outward wave by the artificial boundary.The theoretical derivation of the wave function combination method for the 3D problem and the validation of the program are introduced，as are the basic steps needed to apply this method to the earthquake input.

high concrete-faced rockfill dam； three dimension； non-uniform； displacement input

2014-08-20

國(guó)家自然科學(xué)基金項(xiàng)目(51038007,51079074,51129902)

姚 虞(1989-)，男，博士生，研究課題為高壩抗震計(jì)算.E-mail：yaoyutsinghua@126.com

P315.31

A

1000-0844(2015)02-0324-005

10.3969/j.issn.1000-0844.2015.02.0324