地下結(jié)構(gòu)地震響應(yīng)的計算模型1)

林 皋

(大連理工大學(xué)海岸與近海工程國家重點實驗室，大連116024)(大連理工大學(xué)工程抗震研究所，大連116024)

地下結(jié)構(gòu)地震響應(yīng)的計算模型1)

林 皋2)

(大連理工大學(xué)海岸與近海工程國家重點實驗室，大連116024)(大連理工大學(xué)工程抗震研究所，大連116024)

地震時地下結(jié)構(gòu)在圍巖的約束作用下發(fā)生變形，其動態(tài)特性與地面結(jié)構(gòu)有很大不同.自二十世紀(jì)七八十年代以來，地下結(jié)構(gòu)抗震設(shè)計與研究取得了很大進展.總的看來，工程設(shè)計中普遍采用的計算方法與設(shè)計導(dǎo)則大都建立在比較簡單假定的基礎(chǔ)上，實際的巖土介質(zhì)條件都是十分復(fù)雜的.地下結(jié)構(gòu)抗震研究的近期成果則表現(xiàn)在對地下結(jié)構(gòu)動力分析中的波動散射問題提出了波函數(shù)展開法以及邊界積分方程方法等多種計算方法.但計算相對復(fù)雜，在工程設(shè)計中的推廣應(yīng)用有一定困難.此文致力于地下結(jié)構(gòu)計算模型的改進，使之具有良好的計算精度與效率，又便于工程應(yīng)用.為此，提出了一種地下結(jié)構(gòu)抗震響應(yīng)分析的新的計算模型.模型具有較廣泛的適應(yīng)性，可以進行河谷、孔洞、地下鐵道、隧洞等地下結(jié)構(gòu)的散射與繞射分析.對于復(fù)雜層狀的地質(zhì)條件，提出了格林函數(shù)求解簡便而有效的方法.數(shù)值算例論證了方法的精度和效率.

地下結(jié)構(gòu)，地震響應(yīng)，波的散射，格林函數(shù)，層狀半空間地基

引言

地下工程結(jié)構(gòu)在城市生產(chǎn)和生活中發(fā)揮著重要作用.地下結(jié)構(gòu)承擔(dān)著城市公共設(shè)施的許多重要功能，在通訊、能源供應(yīng)(供電、供氣)、供水排水、交通運輸(地下鐵道、海底隧道)、國防和民防(防空、防爆)等許多方面，發(fā)揮著現(xiàn)代工業(yè)生產(chǎn)和城鎮(zhèn)生活中大動脈的作用.地下結(jié)構(gòu)還常用作儲油庫、儲氣罐；地下結(jié)構(gòu)也較常用來建設(shè)地下電站的廠房.隨著國民經(jīng)濟和城市建設(shè)的發(fā)展，地下工程結(jié)構(gòu)的高度和跨距也愈來愈大，這對地下結(jié)構(gòu)的抗震設(shè)防提出了更高的要求.

歷史上多次大地震中地下結(jié)構(gòu)的破壞造成城市供電、供水、通訊的中斷，引起國民經(jīng)濟的重大損失.著名事例[1]如：1906年美國舊金山大地震(震級 8.2)，由于供水中斷造成火災(zāi)蔓延，破壞區(qū)域達12.2km2(4.7平方哩)，大火燃燒數(shù)日，所造成的損失占地震總損失的80%；1923年日本關(guān)東大地震(震級8.2)，在東京80%的城區(qū)中引起火災(zāi)，并有25條鐵路隧洞受有破壞，橫濱市由于排水總管破壞，除火災(zāi)外又引起次生洪水災(zāi)害.由于供電、供水等地下管線的破壞對民眾的健康和安全關(guān)系重大，因此，這些地下工程結(jié)構(gòu)又稱為生命線工程.

20世紀(jì)80年代以前，國際上還缺乏地下結(jié)構(gòu)或生命線工程抗擊地震荷載的設(shè)計規(guī)范化條款.對一些重要的地下管道，日本采用震度法(地震系數(shù)法)進行設(shè)計，這并不能反映地下結(jié)構(gòu)所遭受地震作用的實際情況.美國舊金山大地震(1906)、日本關(guān)東大地震 (1923)、新潟地震 (1964)、宮城縣地震(1978)和日本海中部地震(1983)中管道等地下結(jié)構(gòu)遭到很大程度的破壞，給城市的生產(chǎn)和生活造成了很大影響，促進了地下結(jié)構(gòu)抗震研究的發(fā)展.這一時期(主要是二十世紀(jì)七八十年代)發(fā)表了大量研究論文[2]，涵蓋地下結(jié)構(gòu)地震震害和地震響應(yīng)觀測、影響因素、地下結(jié)構(gòu)地震反應(yīng)分析和計算模型的建立、設(shè)計要領(lǐng)和設(shè)計準(zhǔn)則等各個方面.與此同時，日本進行了地下結(jié)構(gòu)地震響應(yīng)的觀測和模型試驗研究，其中比較有代表性的是田村和久保等的工作[37].

1977年 8月美國土木工程師學(xué)會生命線地震工程學(xué)科技委員會召開了地震工程學(xué)專業(yè)會議[2]，1981年第7屆世界地震工程會議期間舉行了生命線地震工程學(xué)最新進展委員會的會議.對地下結(jié)構(gòu)或生命線工程抗震研究的情況進行了總結(jié)分析，深化了地下結(jié)構(gòu)抗震特性的認(rèn)識.

根據(jù)會議總結(jié)和會議相關(guān)論文的論述[27]，地下結(jié)構(gòu)的地震響應(yīng)具有以下特點：

(1)大多數(shù)現(xiàn)場觀測資料表明，地震時地下管線和水下隧道無論軸向和側(cè)向均隨圍巖發(fā)生著相同的運動.

(2)觀測到的地下管線和水下隧道的軸向應(yīng)變遠較彎曲應(yīng)變顯著.管道轉(zhuǎn)彎處的彎曲變形與直線段處的彎曲變形具有相同的量級.隧道的軸向剛度相對較大，對圍巖變形起一定的限制作用.

(3)地下管線和隧道地震時產(chǎn)生的慣性力，對結(jié)構(gòu)自身的地震響應(yīng)只產(chǎn)生非常小的影響.

(4)地下管線的存在對圍巖地震動的特性和擾動只有微弱的影響.

本文只討論由地震波傳播所引起的地下結(jié)構(gòu)的地震響應(yīng)分析.

1 地下結(jié)構(gòu)地震響應(yīng)計算模型的發(fā)展

由于地震時地下結(jié)構(gòu)受圍巖介質(zhì)的約束作用而發(fā)生變形，引起震害，因此地下結(jié)構(gòu)的地震響應(yīng)分析應(yīng)包括以下三方面的內(nèi)容：

(1)圍巖介質(zhì)，特別是復(fù)雜巖土介質(zhì)中自由場的地震波動特性.由于地下結(jié)構(gòu)的存在對總體地震波動場的擾動一般較小，因此，可首先研究地下結(jié)構(gòu)建造以前自由場的波動特性.

(2)地震時地下結(jié)構(gòu)周圍的地震波場特性.包括自由場和散射場兩部分.重點需要計算的是由于地下結(jié)構(gòu)不同的幾何形狀和動力剛度特性引起的散射波場特性.

(3)地震時地下結(jié)構(gòu)和周圍巖土介質(zhì)間的動力相互作用特性.由于地下結(jié)構(gòu)的質(zhì)量和剛度與周圍巖土介質(zhì)相比是一個微量，因此地下結(jié)構(gòu)自身的慣性力對其地震變形只起很微小的作用.但是地下結(jié)構(gòu)的剛度，尤其是軸向剛度仍然對圍巖的地震變形發(fā)生足夠的反作用，這使地震時地下結(jié)構(gòu)和圍巖介質(zhì)間產(chǎn)生有差異的地震變形.

二十世紀(jì)七八十年代以來，地下結(jié)構(gòu)的抗震分析取得了很大的進展，但是，對上述三方面的研究深度，仍然不能使人滿意.下文將對計算模型、規(guī)范標(biāo)準(zhǔn)以及文獻中較廣泛采用的有限元數(shù)值計算方法等方面，根據(jù)筆者的體會，作簡要闡述.

自1972年Trifunac[8]發(fā)表關(guān)于半圓形河谷在SH平面波作用下的散射效應(yīng)的論文以來，文獻中關(guān)于地下結(jié)構(gòu)地震響應(yīng)的計算模型方面發(fā)展了很多數(shù)值計算方法.主要針對不同形狀的河谷、沖積層、地下孔洞、地下管道、地下隧道等在SH波、SV波、P波和瑞利波作用下的散射效應(yīng)提出了許多計算模型與方法.具體可參見有關(guān)文獻闡述[912].比較有代表性并獲得較廣泛應(yīng)用的數(shù)值計算方法主要有波函數(shù)展開法和邊界積分方程方法兩種.

波函數(shù)展開(WFE)法主要適用于均質(zhì)空間或半空間散射問題的求解[1316].大多數(shù)求解的為二維散射問題，三維問題的求解局限于軸對稱散射[17]，個別的也求解了非軸對稱散射問題[18].

邊界積分方程(BIE)方法又稱邊界元方法，其特點是只需進行散射體表面的離散，使問題降階一維；同時滿足無窮遠處的輻射條件.但需獲得計算域的基本解，增加了問題求解的困難，并增加了計算工作量.邊界積分方程方法又可區(qū)分為直接BIE方法和間接BIE方法兩種[1925].

將邊界積分方程方法與Haskell和Thomson[26-27]提出的傳遞矩陣方法或Kausel[28]提出的剛度矩陣方法相結(jié)合，可以求解層狀地基的格林函數(shù)，從而求解層狀半空間的散射問題.

文獻中也有的將邊界積分方程法與格林函數(shù)的離散波函數(shù)展開法相結(jié)合進行散射問題的求解[2930].一般情況下，邊界積分方程不適于求解復(fù)雜不均勻介質(zhì)的散射問題，而有限元法則具有較廣泛的適應(yīng)性，對計算域介質(zhì)的不均質(zhì)問題處理非常方便，將兩者相結(jié)合可以發(fā)揮各自的優(yōu)越性.近場采用有限元法分析，遠場輻射邊界則采用邊界積分方程求解[31].同理，也可以將有限元法與波函數(shù)展開法結(jié)合求解散射問題[32].

近年來，中國學(xué)者將前蘇聯(lián)學(xué)者Muskhelishvili所提出的求解彈性力學(xué)問題的復(fù)變函數(shù)法[33]應(yīng)用于求解散射問題也取得了一定效果[3435].此外，文獻中還提出了一些其他方法，本文將不再贅述.

二十世紀(jì)七八十年代以后，一些國家的抗震設(shè)防標(biāo)準(zhǔn)，主要是土木工程的抗震規(guī)范和核電站的抗震規(guī)范,包含了地下結(jié)構(gòu)的抗震設(shè)計準(zhǔn)則.但相對于房屋和橋梁等地面結(jié)構(gòu)來說，條文較簡略，難以滿足實際抗震設(shè)計的要求[3638].

現(xiàn)有規(guī)范中關(guān)于地下結(jié)構(gòu)的抗震分析方法，主要針對地下埋設(shè)管道的計算.對于無限長管，假設(shè)管道的地震變形與均勻無限介質(zhì)中波傳播產(chǎn)生的變形完全相同[3940]，據(jù)此計算管的軸向應(yīng)變與彎曲應(yīng)變.這一做法在核電結(jié)構(gòu)等規(guī)范[41]中應(yīng)用較為普遍.

一般規(guī)范標(biāo)準(zhǔn)中，很少有考慮復(fù)雜地基介質(zhì)中的波傳播特性對地下結(jié)構(gòu)地震響應(yīng)影響的規(guī)定.只有日本沉埋隧道的抗震設(shè)計規(guī)程提出了考慮沿線地基地質(zhì)條件變化的地震動輸入計算模型[38].其中用多質(zhì)點彈簧阻尼體系來模擬圍巖介質(zhì)中地震波的傳播.

文獻中應(yīng)用有限元方法進行地下結(jié)構(gòu)的抗震分析較普遍.將地下結(jié)構(gòu)與圍巖介質(zhì)劃定一范圍作為計算域，用有限元離散進行分析.這種方法可以考慮地下結(jié)構(gòu)的各種幾何形狀和圍巖介質(zhì)的不同特性.但是由于缺乏有關(guān)地下結(jié)構(gòu)抗震設(shè)計的導(dǎo)則和標(biāo)準(zhǔn)，因此，在網(wǎng)格的劃分和布置方面比較自由，不同研究者的計算結(jié)果離散性較大.此外，在計算域范圍、邊界條件的選擇等方面都含有一定的任意性.例如，側(cè)面邊界有的采用能量傳遞邊界，有的采用黏性邊界，有的則采用水平向可自由滑行的邊界；底面邊界有的采用黏性邊界，有的則采用一定深度處的固定邊界.地震動輸入的模型也有一定的任意性.總體上，對自由場地震動的分析，散射場的計算以及地下結(jié)構(gòu)與圍巖介質(zhì)的動力相互作用分析方面是否與實際情況相符都缺乏必要的檢驗，計算精度及其可靠性缺乏論證.將這樣的計算結(jié)果用來進行地下結(jié)構(gòu)的抗震設(shè)計與抗震安全評價是難以令人滿意的.

地下空間的利用發(fā)展很快，適應(yīng)各種需要的地下結(jié)構(gòu)的高度、跨距等都逐漸增加，這對地下結(jié)構(gòu)的抗震設(shè)計提出了新的要求.提高地下結(jié)構(gòu)抗震研究的水平勢在必行.現(xiàn)有的地下結(jié)構(gòu)地震觀測的資料還比較缺乏，并具有一定的局限性，多偏重于可能發(fā)生的地震破損形態(tài)，例如，沿隧道全線或斷面內(nèi)若干測點的變形和加速度等，對地震作用下不同地下結(jié)構(gòu)的變形特點尚缺乏深入認(rèn)識.用以進行新形勢下地下結(jié)構(gòu)的地震響應(yīng)分析還有待進一步深入研究.

2 建議的地下結(jié)構(gòu)地震響應(yīng)分析的基本計算模型

本文的目的在于尋求一種理論上容易理解，計算不復(fù)雜，具有較好的計算精度，同時又便于在實際工程中應(yīng)用的新方法，以便為地下結(jié)構(gòu)抗震設(shè)計導(dǎo)則和規(guī)程的制訂提供一定的參考.

下面從土與結(jié)構(gòu)動力相互作用 (dynamic soilstructure interaction)的基本理論出發(fā)來建立地下結(jié)構(gòu)地震響應(yīng)分析的方程[42].應(yīng)用子結(jié)構(gòu)分析方法，土與結(jié)構(gòu)動力相互作用的計算方程可表示如下(圖1)

式中,ω為激勵頻率,ζ表示滯回阻尼系數(shù).如需考慮高階阻尼影響時，S的表達式可參見文獻[43],進行地震作用分析時Ps=0，結(jié)構(gòu)和無限地基的相互作用力Rb可表示為

圖1 結(jié)構(gòu)與無限地基系統(tǒng)Fig.1 Structure-unbounded soil system

將式(2)和式(3)代入式(1)可得

在圖 1(b)中，將地基被挖去的部分表示為結(jié)構(gòu)，則式(5)中矩陣Sbb-SbsSSsb表示地基被挖去部分的剛度(指被挖去部分凝聚于邊界b上的剛度).這樣式(5)左方括號內(nèi)的矩陣+即表示地基未開挖前b邊界的動力剛度；而即表示地基未開挖前b邊界的自由場地震動[44].

圖2 散射體的計算模型示意圖Fig.2 Computational model for surface and subsurface irregularities and inhomogeneities

需要指出，式(8)的應(yīng)用需要根據(jù)情況作出相應(yīng)的變化.例如對于沖擊層填充的河谷(圖3)，計算公式相應(yīng)地轉(zhuǎn)換為

圖3 填充河谷散射體Fig.3 Scattering by alluvial valley

圖4 散射場計算邊界的選擇Fig.4 Selection of the outer boundaries for scattering analysis

可以指出，這里所建立的計算公式，既可適用于二維結(jié)構(gòu)，也可適用于三維結(jié)構(gòu)；既可適用于均質(zhì)半

3 數(shù)值驗算

本文將以均質(zhì)半無限空間中的散射問題為主，通過算例檢驗計算模型的有效性.為計算均質(zhì)半無限空間的動力剛度,比例邊界有限元方法(SBFEM)[45]是有效的方法.SBFEM為半解析半數(shù)值性的求解方法，其徑向可獲得解析解，環(huán)向則達到有限元的計算精度.由于只需要進行邊界面的離散，問題的維數(shù)降低一階；并可自動滿足無限遠處的輻射條件.因此SBFEM兼有FEM和BEM的優(yōu)點，也不需要基本解，計算簡便.

3.1 半圓形河谷的散射

如圖5所示,河谷邊界ABC,計算邊界DEFG，相應(yīng)計算域大小為4a×2a(a為河谷半徑).設(shè)圍巖介質(zhì)材料特性參數(shù)為：質(zhì)量密度ρ0=1.5，剪切波速cs0=1，泊松比ν0=1/3；河谷中填充材料參數(shù)為：質(zhì)量密度ρ1=1，剪切波速cs1=0.5，泊松比ν0=1/3.無量綱頻率 η =(ωa)/(πcs)，波長 λ =2/η.計算采用式(9).網(wǎng)格離散如圖5.各子域邊界采用了3節(jié)點二次單元離散，共144單元，288節(jié)點.相應(yīng)DEFG邊界設(shè)48單元，97節(jié)點.入射波從底部垂直輸入，計算結(jié)果精度高，與Luco和de Barros的解[46]相符性良好 (圖 6).Luco和 de Barros，以及 Dravinski和Mossessian[47]采用的都是間接邊界積分方程(IBIE)方法.

圖5 半圓形河谷的計算圖形Fig.5 Schematic view of semi-circular valley

圖6 半圓形沉積河谷散射問題的檢驗Fig.6 Verificatio of the scattering by semi-circular valley

圖6 半圓形沉積河谷散射問題的檢驗(續(xù))Fig.6 Verificatio of the scattering by semi-circular valley(continued)

3.2 地下洞室襯砌的地震響應(yīng)

研究地下隧道襯砌的地震響應(yīng).de Barros和Luco[48]求解了這一問題，采用間接邊界積分方程方法模擬外圍半空間土介質(zhì)的作用，再結(jié)合Donnell的殼體理論計算隧洞襯砌外半徑的波動響應(yīng).子域括號內(nèi)的數(shù)字代表該子域相似中心的坐標(biāo).隧道結(jié)構(gòu)如圖7(a)所示.隧道埋深H=8.33r1=7.573a.襯砌材料質(zhì)量密度 ρ0=2.24×103kg/m3，彈性模量E0=1.6×1010N/m2，泊松比 ν0=0.2.圍巖地基介質(zhì)材料質(zhì)量密度ρ1=7.665×103kg/m3，彈性模量 E1=6.9×108N/m2，泊松比 ν0=0.45.位移按r=a進行規(guī)格化.假設(shè)SV波垂直入射，無量綱頻率η=ωa/(πcs)=0.132.計算的襯砌響應(yīng)如圖8所示.可見與文獻解的相符性較好.

圖7 均質(zhì)半無限地基中的隧道襯砌Fig.7 Tunnel and its liner embedded in uniform half-space

圖8 SV波垂直入射作用下隧道襯砌外表面的動力響應(yīng)Fig.8 Dynamic response at the outer radius of the tunnel liner subjected to vertically incident SV wave

4 復(fù)雜地基條件下地下結(jié)構(gòu)的地震響應(yīng)分析

其中的水平分層地基比較常見，下文將重點研究.對于分塊不均質(zhì)地基，當(dāng)各塊之間可以用同一相似中心加以描繪時(圖9(b))，提出采用比例邊界有限元的處理方法[49]，對近場含不規(guī)則不均質(zhì)捕虜體等復(fù)雜地基(圖9(c))則提出采用分步阻尼影響抽取法[50]進行求解.所提出的方法雖然是結(jié)合結(jié)構(gòu)與地基動力相互作用分析進行闡述，但可適用于地下結(jié)構(gòu)的抗震分析.

圖9 不均質(zhì)地基的計算模型Fig.9 Computational model for inhomogeneous unbounded soil

對各向同性和非各向同性層狀半空間不均質(zhì)地基的動力剛度及格林函數(shù)的求解，作者研究組從2012年起提出了基于積分變換的方法，并逐步得到完善[43,5153]，并可適用于橫觀各向同性的層狀介質(zhì)地基的求解.

下面介紹其基本思想和主要方程.

4.1 層狀介質(zhì)的波動方程

在圓柱坐標(biāo)系下，層狀介質(zhì)(圖10)中任一層以位移所表示的波動方程具有如下形式

圖10 水平層狀半空間Fig.10 Horizontally layered strata overlying half-space

對式(10)進行Fourier-Bessel變換，得出波數(shù)域的表達式.變換式如下

相應(yīng)的反變換式為

式中k代表波數(shù)，其余系數(shù)的表達式為

式(15)中的兩個子式，上面的子式對應(yīng)于荷載與變形相對于X軸對稱的情況；下面的子式則為荷載與變形相對于X軸成反對稱的情況.其中Jn(kr)為柱貝塞爾函數(shù).Green函數(shù)計算時，針對圓形單元面上作用的均布水平向荷載和豎向荷載，可選擇n=1和n=0兩種情況(圖11)與之對應(yīng).當(dāng)圓形單元半徑趨于極限Δr→0時，均布荷載轉(zhuǎn)化為集中水平力和集中豎向力作用的情況.即，在實際應(yīng)用時，主要可限于n=1和n=0的兩種情況.

圖11 Green函數(shù)計算荷載單元Fig.11 Loaded element for evaluation of Green’s function

按式(12)進行Fourier-Bessel變換后得出頻率--波數(shù)域內(nèi)解耦的SV-P波動方程(變量為ur.uz)和SH波動方程(變量為uθ)[5152]，也可參見文獻[28].

對橫觀各向同性介質(zhì)中的波動，解耦方程同樣成立，但相應(yīng)系數(shù)改變?nèi)缦?/p>

式中,上標(biāo)m=1代表SV-P波動，m=2代表SH波動.Im為單位陣，m=1為2階單位陣，m=2為1階單位陣.

對各向同性介質(zhì)(方程(17)與方程(18))

對于橫觀各向同性介質(zhì)

為了進行微分方程(21)的求解，引入如下定義的應(yīng)力矢量pm作為位移um的對偶變量

可以證明pm滿足關(guān)系式[54]

以下的求解過程對m=1和m=2都適用，所以在下文中可將m省略.

利用式(26)，微分方程(21)可化為如下對耦形式

矩陣A,B,C,D和矩陣Kij(i,j=1,2)的變換關(guān)系如下

在狀態(tài)空間，對耦微分方程(27)可以聯(lián)立求解

齊次一階線性常微分方程(29)的解為指數(shù)函數(shù).

采用鐘萬勰提出的精細積分方法[55]可以將解進行精確計算.

對層狀半空間中的任意一層，或其中的一層，設(shè)層厚為η=zb-za，則層上、下兩端的位移和應(yīng)力將成立如下關(guān)系

或?qū)憺?/p>

式中

為達到高精度的計算效果，T按如下方式進行計算[55]

式中，τ=η/b，b可取為任意整數(shù).當(dāng)取b=2N,N=20時，這相當(dāng)于將η分成2N=1048557個微細薄層用式(34)計算，從而可以達到任意希望的精度，而計算工作量不大.文獻[55]中給出了計算T的遞推公式，應(yīng)用方便.

為了進行層的結(jié)合和計算Green函數(shù)，解也宜表示成如下對耦變量的形式

式中矩陣F,G,Q也由T求出

對于水平分層介質(zhì)，采用積分變換方法，可以求得波動方程在波數(shù)域的簡便形式(式(17)和式(18)或式(19)和式(20)).進一步引入對偶變量，又可將其化為一階齊次線性常微分方程(29)，其解可采用精細積分方法達到任意希望的精度.最后將解表示成對耦形式(36)，可便于分層之間的結(jié)合和Green函數(shù)的求解，說明如下.

4.2 相鄰層的結(jié)合

當(dāng)將相鄰的兩層：層1[za,zb]和層2[zb,zc]進行結(jié)合時，應(yīng)用式(36)，可有

從式(38)中消去ub和pb，即可求得合并后新層c[za,zc]的力與位移關(guān)系

式中

4.3 Green函數(shù)的計算

層狀半空間內(nèi)部點Green函數(shù)的計算是地下結(jié)構(gòu)抗震分析的核心部分.文獻中廣泛采用的剛度矩陣方法(sti ff ness matrix method)需要進行大型矩陣求解[28].這里提出對偶變量的轉(zhuǎn)換方法，計算簡捷方便，同時又準(zhǔn)確高效.計算圖形如圖12和圖13所示.圖中源點(source)表示力的作用點，接受點(receiver)表示位移的計算點.已知源點和接受點的位置，可將層狀半空間劃分為3個子集，應(yīng)用4.2節(jié)中闡述的方法可以分別建立各個子集形如式(36)所示的對耦變量方程，并相應(yīng)求出其Fi,Qi,Gi矩陣(i=u,m,l，其中u表示上部，m表示中部，l表示下部)方程.Green函數(shù)的計算，基本上可劃分為圖13所示的6種工況.應(yīng)用式(36)，即可求得各種工況的Green函數(shù)計算的矩陣表達式[53].

圖12 層狀半空間的Green函數(shù)Fig.12 Green’s function for multilayered half-space

圖13 層狀半空間不同工況的Green函數(shù)計算Fig.13 Computation of Green’s function for various cases of boundary conditions of multilayered half-space

下面給出不同工況的Green函數(shù)的矩陣表達式.下標(biāo)a,b表示源點和接受點所在位置.

(1)工況(a)和(b)層狀地基表面點的Green函數(shù)

工況(a)下部為剛性地基

工況(b)下部為彈性半空間

(2)工況(c)和(d)下部固定的層狀地基內(nèi)部點的格林函數(shù)

工況(c)

(3)工況(e)和(f)下部為彈性半空間的層狀地基內(nèi)部點的格林函數(shù)

工況(e)

所有矩陣的階數(shù)對SV-P波為2×2，對SH波為1×1，計算很簡便.最后可將波數(shù)域各種工況格林函數(shù)的計算，寫成如下形式的統(tǒng)一的表達式(參見式(17)和式 (18)或式(19)和式(20),ur,uz和uB是解耦的，可分別求解

對圖 11表示的格林函數(shù)計算的荷載形式，例如，半徑為Δr的圓形單元上分布的Z向和X向均布荷載和，其相應(yīng)的波數(shù)域的表達式為

在上述式中，令 Δr→0，πΔr2z→z或πΔr2x→x就可得出集中力z或x作用的表達式.再將式(53)、式(54)代入式(51)、式(52)并利用(13)進行反變換，就可得出圓柱坐標(biāo)系物理域中格林函數(shù)的表達式，可統(tǒng)一寫成

4.4 數(shù)值檢驗

研究均質(zhì)半無限地基上表面彈性層中半圓形河谷的散射(圖14).散射計算需要求解層狀半空間中的格林函數(shù).計算中介質(zhì)材料特性參數(shù)假設(shè)為：均質(zhì)半空間 ρ0=4/3，r0=0.333，cs0=2；表面彈性層 ρ1=1，r1=0.333，cs1=1；河谷中沖積層ρ2=0.667，r2=0.333，cs2=0.5.幾何尺寸：河谷半徑r，表面彈性層厚h=1.5r.無量綱頻率η=ωa/(πcs).

在SV波垂直入射情況下半圓形充填河谷地表位移幅值計算結(jié)果如圖15所示，可見與文獻[56]的相符性很好.

圖14 半無限地基上表面彈性層中的半圓形河谷的散射Fig.14 Scattering of a semi-circular valley embedded in the surface layer overlying elastic half-space

圖15 SV波垂直入射情況下層狀半空間中半圓形河谷地表位移幅值Fig.15 Amplitude of displacement response of the semi-circular valley embedded in layered half-space under vertically incident SV wave

5 地下結(jié)構(gòu)地震響應(yīng)分析的其他有關(guān)問題

5.1地下結(jié)構(gòu)地震響應(yīng)的時域分析

地下結(jié)構(gòu)地震響應(yīng)的時域計算應(yīng)用運動方程(1)，關(guān)鍵是相互作用力 R=或R=的計算.由于無限域中R是激勵頻率的函數(shù)，需將其轉(zhuǎn)換成時域的表達式，如式(51)所示(假設(shè)初值條件

式中，S∞(ω)代表無限域的動剛度(ω)或ω)；S∞(t)代表單位位移脈沖響應(yīng)函數(shù)

式(58)只具有形式上的意義.由于S∞(ω)非平方可積，含有奇異性，式(58)不能直接用于脈沖響應(yīng)函數(shù)的計算[57]，需進行必要的處理，將另文闡述.

5.2 地下結(jié)構(gòu)地震響應(yīng)的簡化計算

地下結(jié)構(gòu)也可采用簡化的數(shù)值方法進行初步的抗震分析.以圖16所示的地下結(jié)構(gòu)為例來加以說明.圖中虛線范圍以內(nèi)為計算域，域內(nèi)可采用有限元等數(shù)值方法進行分析.根據(jù)第2節(jié)計算模型的推導(dǎo)，虛線所表示的邊界應(yīng)滿足相互作用力的條件R=或R=.也可選擇較大的計算域，并設(shè)立能量傳遞邊界，通過邊界進行自由場的地震動輸入.文獻中一般的能量傳遞邊界主要建立在均質(zhì)無限地基假定的基礎(chǔ)上，不適于復(fù)雜地基條件下地下結(jié)構(gòu)的地震響應(yīng)分析.如采用時間、空間非耦合的局部人工透射邊界，計算準(zhǔn)確度并不理想，需要采用高階局部人工邊界，但高階局部邊界計算精度的提高是以計算自由度的增加為代價的.我們提出的層狀地基耦合性(或全域性)時域傳遞邊界[58]，不僅計算精度高，而且計算效率也高，可供參考.

圖16 地下結(jié)構(gòu)抗震分析的簡化模型Fig.16 A simplifie model for seismic analysis of under-ground structures

6 結(jié)束語

提出一種基于結(jié)構(gòu)與地基動力相互作用方程的地下結(jié)構(gòu)抗震分析模型，可以很方便地進行地震波散射場與地下結(jié)構(gòu)和圍巖動力相互作用的計算.模型具有廣泛的適用性.可以考慮層狀半空間等復(fù)雜地基的影響.計算準(zhǔn)確、高效.本文以二維均勻半空間中地下結(jié)構(gòu)的地震響應(yīng)為例，闡述地下結(jié)構(gòu)抗震計算模型的設(shè)想.但是對復(fù)雜層狀半空間地下結(jié)構(gòu)的地震響應(yīng)只提出了計算方法和相關(guān)方程還來不及展開，以后將陸續(xù)進行闡述.

致謝 博士生李志遠、韓澤軍和課題組李建波、胡志強等在本文完成中發(fā)揮了重要作用,在此表示感謝.

1 Ariman T,Gregory E Muleski.A review of the response of buried pipelines under seismic excitations.Earthquake Engineering&Structural Dynamics,1981,9(2)：133-152

2 Wang LRL,O’Rourke MJ.Overview of buried pipelines under seismic loading.Journal of Technical Councils of ASCE,1978,104(TC1)：121-130

3 Okamoto S,Tamura C.Behavior of subaqueous tunnel during earthquakes.International Journal of Earthquake Engineering&Structural Dynamics,1973,1(3)：253-266

4 Tamura C,Okamoto S,Kato K.Observations on dynamic strains of submerged tunnel during earthquakes.Bulletin of Earthquake Resistant Structure Research Center,University of Tokyo,1972,6：1-21

5 Tamura C.Dynamic behavior of submerged tunnel during earthquake.Transactions of the Japan Society of Civil Engineers,1975,7：197-200

6 Kubo K.Behavior of underground water pipes during an earthquake//Proceedings of the Fifth World Conference on Earthquake Engineering,Rome,Italy,1974,569-578

7 Kubo K,Katayama T,Ohashi A.Present state of lifeline earthquake engineering in Japan//The Current State of Knowledge of Lifeline Earthquake Engineering,Rome,Italy,1974,569-578

8 Trifunac MD.Scattering of plane SH waves by a semi-cylindrical canyon.Earthquake Engineering&Structural Dynamics,1972,1(3)：267-281

9 Luco JE,De Barros FCP.Seismic response of a cylindrical shell embedded in a layered viscoelastic half-space.I：Formulation.Earthquake Engineering&Structural Dynamics,1994,23(5)：553-567

10 Wong HL.E ff ects of surface topography and site conditions//Proc.of the EPRI Workshop on Strong Ground Motion Simulation and Earthquake Engineering Applications,Los Altos,California,April 30-May 3,1984

11 Sanchez-Sesma FJ.Site e ff ects on strong ground motion.Soil Dynamics and Earthquake Engineering,1987,6(2)：124-132

12 Von Thun JL.Earthquake Engineering and Soil Dynamics II—Recent Advances in Ground-Motion Evaluation.[s.l.]：ASCE,1988

13 El-Akily N,Datta SK.Response of a circular cylindrical shell to disturbances in a half-space.Earthquake Engineering&Structural Dynamics,1980,8(5)：469-477

14 Wong KC,Shah AH,Datta SK.Dynamic stresses and displacements in a buried tunnel.Journal of Engineering Mechanics,1985,111(2)：218-234

15 Eshraghi H,Dravinski M.Transient scattering of elastic waves by dipping layers of arbitrary shape part 1：Antiplane strain model.Earthquake Engineering&Structural Dynamics,1989,18(3)：397-415

16 Eshraghi H,Dravinski M.Transient scattering of elastic waves by dipping layers of arbitrary shape.Part 2：Plane strain model.Earthquake Engineering&Structural Dynamics,1989,18(3)：417-434

17 Lee VW.Three-dimensional di ff raction of plane P,SV&SH waves by a hemispherical alluvial valley.International Journal of Soil Dynamics and Earthquake Engineering,1984,3(3)：133-144

18 Eshraghi H,Dravinski M.Scattering of plane harmonic SH,SV,P and Rayleigh waves by non-axisymmetric three-dimensional canyons：A wave function expansion approach.Earthquake Engineering&Structural Dynamics,1989,18(7)：983-998

19 Wong HL,Jennings PC.E ff ects of canyon topography on strong ground motion.Bulletin of the Seismological Society of America,1975,65(5)：1239-1257

20 Zhang L,Chopra AK.Three-dimensional analysis of spatially varying ground motions around a uniform canyon in a homogeneous half-space.Earthquake Engineering&Structural Dynamics,1991,20(10)：911-926

21 Wong HL.E ff ect of surface topography on the di ff raction of P,SV,and Rayleigh waves.Bulletin of the Seismological Society of America,1982,72(4)：1167-1183

22 S′anchez-Sesma FJ,Campillo M.Di ff raction of P,SV,and Rayleigh waves by topographic features：A boundary integral formulation.Bulletin of the Seismological Society of America,1991,81(6)：2234-2253

23 Pedersen HA,Sanchez-Sesma FJ,Campillo M.Three-dimensional scattering by two-dimensional topographies.Bulletin of the Seismological Society of America,1994,84(4)：1169-1183

24 Luco JE,Wong HL,De Barros FCP.Three-dimensional response of a cylindrical canyon in a layered half-space.Earthquake Engineering&Structural Dynamics,1990,19(6)：799-817

25 Luco JE,De Barros FCP.Seismic response of a cylindrical shell embedded in a layered viscoelastic half-space.I：Formulation.Earthquake Engineering&Structural Dynamics,1994,23(5)：553-567

26 Haskell NA.The dispersion of surface waves on multilayered media.Bulletin of the Seismological Society of America,1953,43(1)：86-103

27 Thomson WT.Transmission of elastic waves through a stratifie solid medium.Journal of Applied Physics,1950,21(2)：89-93

28 Kausel E.Fundamental Solutions in Elastodynamics：A Compendium.New York,USA：Cambridge University Press,2006

29 Bouchon M.A simple,complete numerical solution to the problem of di ff raction of SH waves by an irregular surface.The Journal of the Acoustical Society of America,1985,77(1)：1-5

30 Kawase H.Time-domain response of a semi-circular canyon for incident SV,P,and Rayleigh waves calculated by the discrete wavenumber boundary element method.Bulletin of the Seismological Society of America,1988,78(4)：1415-1437

31 Mossessian TK,Dravinski M.Application of a hybrid method for scattering of P,SV,and Rayleigh waves by near-surface irregularities.Bulletin of the Seismological Society of America,1987,77(5)：1784-1803

32 Shah AH,Wong KC,Datta SK.Di ff raction of plane SH waves in a Half-space.Earthquake Engineering&Structural Dynamics,1982,10(4)：519-528

33 Muskhelishvili NI.Some Basic Problems of the Mathematical Theory of Elasticity.4th Edn.Netherlands,Noorho ff,New York,University of Groningen,1963(in Russian)

34 Liu D,Gai B,Tao G.Applications of the method of complex functions to dynamic stress concentrations.Wave Motion,1982,4(3)：293-304

35 Liu Q,Zhao M,Wang L.Scattering of plane P,SV or Rayleigh waves by a shallow lined tunnel in an elastic half space.Soil Dynamics and Earthquake Engineering,2013,49：52-63

36王汝梁.地下管線的震害、抗震驗算、設(shè)計與措施.王珞珈等譯.北京：地震出版社，2007(Wang Ruliang.Buried Lifeline Earthquake Engineering and Related Research.Wang LJ,et al.Trans.Beijing：Seismicity Press,2007(in Chinese))

37侯忠良.地下管線抗震.北京：學(xué)術(shù)書刊出版社，1990(Hou Zhongliang.Earthquake Resistance of Buried Pipelines.Beijing：Academic Press,1990(in Chinese))

38 Earthquake Engineering Committee.The Japan Society of Engineers.Earthquake Resistant Design for Civil Engineering Structures in Japan.1988

39 Sakurai A,Takahashi T.Dynamic stresses of underground pipelines during Earthquake//Proceedings of 4th World Conference on Earthquake Engineering,Chile,1969,811-895

40 Newmark NM.Problems in wave propagation in soil and rock//Proceedings of International Symposium in Wave Propagation and Dynamic Properties of Earth Materials,Albuquerque,New Mexico,1967,7-26

41 Yeh GCK.Seismic analysis of slender buried beams.Bulletin of the Seismological Society of America,1974,64(15)：1551-1562

42 Lin G,Li ZY,Li JB.A substructure replacement technique for the numerical solution of wave scattering problem(to be published)

43 Song C.The scaled boundary finit element method in structural dynamics.International Journal for Numerical Methods in Engineering,2009,77(8)：1139-1171

44 Lin G,Han Z,Li J.Soil–structure interaction analysis on anisotropic stratifie medium.G′eotechnique,2014,64(7)：570-580

45 Wolf JP.Dynamic Soil-Structure Interaction.[s.l.]：Prentice Hall Int.,1985

46 Luco JE,de Barros FCP.Three-dimensional response of a layered cylindrical valley embedded in a layered half-space.Earthquake Engineering&Structural Dynamics,1995,24(1)：109-125

47 Dravinski M,Mossessian TK.Scattering of plane harmonic P,SV,and Rayleigh waves by dipping layers of arbitrary shape.Bulletin of the Seismological Society of America,1987,77(1)：212-235

48 de Barros FCP,Luco JE.Seismic response of a cylindrical shell embedded in a layered viscoelastic half-space.II：Validation and numerical results.Earthquake Engineering&Structural Dynamics,1994,23(5)：569-580

49 Lin G,Du J,Hu Z.Earthquake analysis of arch and gravity dams including the e ff ects of foundation inhomogeneity.Frontiers of Architecture&Civil Engineering in China,2007,1(1)：41-50

50 YinX,LiJ,WuC,etal.ANSYSimplementationofdampingsolvent stepwise extraction method for nonlinear seismic analysis of large 3-D structures.Soil Dynamics&Earthquake Engineering,2013,44(1)：139-152

51林皋，韓澤軍，李建波.層狀地基任意形狀剛性基礎(chǔ)動力響應(yīng)求解.力學(xué)學(xué)報，2012，44(6)：1016-1027(Lin Gao,Han Zejun,Li Jianbo.Solution of the dynamic response of rigid foundation of arbitrary shape on multi-layered soil.Chinese Journal of Theoretical and Applied Mechanics,2012,44(6)：1016-1027(in Chinese))

52 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.Soil Dynamics&Earthquake Engineering,2015,70：48-59

53 Lin G,Han Z,Lu S,et al.Wave motion equation and the Green’s function for transverse isotropic multi-layered half-space.(to be published)

54 Zhong W,Williams FW,Bennett PN.Extension of the wittrickwilliams algorithm to mixed variable systems.Journal of Vibration&Acoustics,1997,119(3)：334-340

55 Zhong WX.On precise integration method.Journal of Computational&Applied Mathematics,2004,163(1)：59-78

56 de Barros FCP,Luco JE.Amplificatio of obliquely incident waves by a cylindrical valley embedded in a layered half-space.Soil Dynamics&Earthquake Engineering,1995,14(3)：163-175

57 Wolf JP,Song C.Finite-element Modelling of Unbounded Media.Chichester：Wiley,1996

58 Lin G,Lu S,Liu J.Transmitting boundary for transient analysis of wavepropagationin layered media formulated based on acceleration unit-impulse response.Soil Dynamics&Earthquake Engineering,2016,90：494-509

A COMPUTATIONAL MODEL FOR SEISMIC RESPONSE ANALYSIS OF UNDERGROUND STAUCTURES1)

Lin Gao2)
(State Key Laboratory of Coastal and O ff shore Engineering，Dalian University of Technology，Dalian 116024，China)(Institute of Earthquake Engineering，Dalian University of Technology，Dalian 116024，China)

Under earthquake excitation,the deformation of underground structures is restricted by the surrounding soil media.The dynamic behavior of it displays quite di ff erently from aboveground structures.Considerable progress has been made on the design and research of the underground structures since seventies-eighties of twentieth century.However,in the main,the widely used computational methods and guidelines in engineering design practice are based on rather simple assumptions,in reality,the actual soil conditions might be much more complex than ideally boundary conditions.Recent achievements of earthquake research on underground structures lie in the development of various computational methods for wave scattering problems of underground structures,such as the wave function expansion method,the boundary integral equation method etc.As the computation is somewhat complex,which impedes its application and dissemination in the engineering design practice.The author devotes himself to the improvement of the computational model for seismic analysis of underground structures,such that it achieves higher accuracy and efficiency,meanwhile it proves to be convenient for engineering design.To this end,a new model for seismic response analysis of underground structuresis proposed.The model is versatile to deal with wave scattering and di ff raction by canyons,subsurface cavities,subways and tunnels etc.In case of the presence of complex soil conditions like the layered half-space,a simple and e ff ective technique is developed for the evaluation of Green’s functions.Numerical examples are provided to validate the accuracy and efficiency of the proposed approach.

underground structure,seismic response,wave scattering,Green’s function,layered half-space.

TU93

：A

10.6052/0459-1879-16-301

2016–10–28 收稿，2017–01–25 錄用,2017–03–02 網(wǎng)絡(luò)版發(fā)表.

1)國家重點研發(fā)計劃資助項目(2016YFB0201001).

2)林皋，中國科學(xué)院院士，教授，主要研究方向：水利工程與核電工程抗震.E-mail：gaolin@dlut.edu.cn

林皋.地下結(jié)構(gòu)地震響應(yīng)的計算模型.力學(xué)學(xué)報,2017,49(3)：528-542

Lin Gao.A computational model for seismic response analysis of underground stauctures.Chinese Journal of Theoretical and Applied Mechanics,2017,49(3)：528-542