三維飽和沉積盆地地震動模擬方法
【專利摘要】本發(fā)明提供一種三維飽和半空間中沉積場地地震動模擬方法，該方法基于間接邊界元方法，推導(dǎo)了全空間飽和斜面圓盤荷載格林函數(shù)，解決了求解邊界的奇異性問題；對于飽和兩相模型模擬固液動力耦合效應(yīng)、孔壓及流量的變化規(guī)律；建立三維飽和場地模型，基于Biot飽和多介質(zhì)理論推導(dǎo)飽和格林函數(shù)，根據(jù)邊界條件建立方程求解，最后進行頻時域地震響應(yīng)模擬分析以期實現(xiàn)未來地震動預(yù)測。本發(fā)明效果是為飽和不均勻介質(zhì)工程波動問題求解提供一種高效算法和實用程序，為沿海復(fù)雜飽和場地區(qū)域地震區(qū)劃、結(jié)構(gòu)抗震設(shè)防提供理論依據(jù)，為實際復(fù)雜工程的模擬帶來顯著的經(jīng)濟效益。
【專利說明】
三維飽和沉積盆地地震動模擬方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明涉及地震工程與工程波動計算領(lǐng)域，特別是一種基于間接邊界元法的三維 飽和半空間中沉積場地地震動模擬方法。
【背景技術(shù)】
[0002] 地震波散射一直是地震工程、地震學(xué)、地球物理等眾多學(xué)科的一個研究熱點。對于 地震波動問題的理論求解方法整體可以分為解析法和數(shù)值法，其中數(shù)值法包括域離散型的 有限元法、有限差分法、離散波數(shù)法和邊界離散的邊界元法等。同域離散型方法相比，邊界 元法(BEM)具有高精度、降低問題求解維度、便于處理高梯度應(yīng)力變化的優(yōu)點，同時特別便 于處理無限域波動問題，無需引入無反射邊界條件，另外還克服了普通有限元法高頻彌散 的難題。邊界元法中由位勢理論建立的間接邊界元法不僅具有高精度的優(yōu)點，且無需經(jīng)過 復(fù)雜數(shù)學(xué)理論推導(dǎo)，易于理解;在編程上易于實現(xiàn)，方法直觀，更適合于任意復(fù)雜問題的求 解。目前，大多采用邊界元法進行研究集中于二維或三維均勻彈性半空間中局部場地對彈 性波的散射問題且已日臻完善。然而，近些年來在我國沿海區(qū)域，超高層建筑、越江隧道、跨 海大橋、以及高速鐵路等大型工程不斷涌現(xiàn)，且很多位于高地震烈度復(fù)雜場地。這些大型工 程的抗震設(shè)計亟需更為精細可靠的設(shè)計地震動參數(shù)。鑒于此，飽和兩相介質(zhì)地震波動模擬 正逐漸成為地震工程及土動力學(xué)領(lǐng)域一個研究前沿。
[0003] 現(xiàn)實中飽和巖土層廣泛存在于地下水位以下、沿海河谷地帶以及壩底淤積等。同 單相固體介質(zhì)相比，由于受固體骨架和孔隙流體的相互作用控制，流體飽和介質(zhì)力學(xué)行為 有其自身特殊性，如慢縱波的存在、孔隙水壓的消長和流量變化，以及地震波驅(qū)動下流體流 動引起衰減及頻散等。通過局部飽和不規(guī)則場地對地震波的散射模擬，以獲得該地不同區(qū) 域的地震動響應(yīng)，最終為城市地震區(qū)域化提供必要的理論依據(jù)。為更好的解釋這些特殊規(guī) 律，提高飽和地層中地震波散射模擬的定量化水平，亟需深入研發(fā)基于飽和兩相介質(zhì)模型 的地震波模擬方法及應(yīng)用研究。
[0004] 鑒于研究問題的復(fù)雜性，目前對于飽和半空間中局部場地對彈性波的散射解答相 對較少，而三維飽和局部場地對彈性波的散射結(jié)果也僅有少數(shù)的解析解。而由于解析法求 解過程的局限性，對于實際工程中的非均勻或較為復(fù)雜的三維場地，解析法將無法精確求 解其對彈性波的散射問題。

【發(fā)明內(nèi)容】

[0005] 本發(fā)明的目的是提供一種三維飽和半空間中沉積場地地震動模擬方法，該方法詳 細推導(dǎo)了全空間飽和斜面圓盤荷載格林函數(shù)，解決了求解邊界的奇異性問題;對于飽和兩 相模型可更好地模擬固液動力耦合效應(yīng)、孔壓及流量的變化規(guī)律。
[0006] 為實現(xiàn)上述目的，本發(fā)明采用的技術(shù)方案是提供一種三維飽和半空間中沉積場地 地震動模擬方法，該方法是在建立三維飽和沉積盆地模型基礎(chǔ)上，通過由飽和土中的土體 骨架和流體之間的相對關(guān)系，推導(dǎo)飽和格林函數(shù)，采用間接邊界元法解決飽和局部場地的 地震動模擬，該方法包括有以下步驟:量化
[0007] (一)建立三維飽和場地模型：
[0008] 考慮平面波P或S和Rayleigh波的入射，基于單層位勢理論，在沉積盆地模型邊界 表面施加三個正交方向的虛擬荷載以構(gòu)造散射波場，采用邊界元方法實施中，需首先將邊 界表面離散為三角形或四邊形單元，為便于處理所述三角形或四邊形單元格林函數(shù)，采用 斜面圓盤均布荷載近似覆蓋在所述邊界表面的離散單元上，即分所述離散單元上虛擬荷載 對自身單元作用為本單元，所述離散單元上虛擬荷載對自身單元之外的離散單元作用為非 本單元，推導(dǎo)出飽和介質(zhì)中三個方向的力F x，F(xiàn)y，F(xiàn)z集中力和流量元Q作用下的格林函數(shù)，進 而根據(jù)所述邊界表面的透水和不透水邊界條件構(gòu)造矩陣方程組：
[0009] [Η] [ Φ ] = [B]，求解得到虛擬荷載密度：[Φ ] = [H^B]
[0010] 其中，[Η]為格林函數(shù)影響矩陣，[Β]為邊界條件解。
[0011] 最后將邊界表面上各個單元上的虛擬荷載的綜合作用而得散射波場和不含局部 不均勻場地時彈性波入射下的半空間自由波場疊加即得到總波場；由于虛擬荷載直接施加 在邊界表面上，能夠解決多種形狀的飽和局部場地。
[0012](二)三維飽和格林函數(shù)推導(dǎo)：
[0013]基于Biot飽和兩相介質(zhì)理論，推導(dǎo)飽和介質(zhì)中三個方向的力Fx，F(xiàn)y，F(xiàn)z集中力和流 量元Q作用下的格林函數(shù)，格林函數(shù)解如下：
[0014] 1非本單元格林函數(shù)解
[0015] 1.1集中力作用
[0016] 類似于三維彈性半空間問題，在均勻各向同性彈性固體介質(zhì)中，位移格林函數(shù)表 達為：
[0017] Gij(x,C) = [f25ij+(fi-f2) γ i γ j]/43iyr (1)
[0018]其中，γ j = (xj-|j)/r，r表示波源點（|1,|2，|3)與接收點（1142 43)之間的距離，即 1'2=(叉1-|1)2+(叉2-|2)2+(叉3-|3)2,51」表示1(1'〇116。1^1(16]^4表示拉梅常數(shù)；1^=〇/0，1^1 =ω /αι，ka2 = ω /(^分別表示S，Pl(快波），p2(慢波)波數(shù);β， αι，(^表示其各自對應(yīng)波速，f 4口 f2定義為：
[0019] ?ι = Αι(β2/αι2) [l-2i/(k0ir)-2/(kair)2]exp(-ikair)+A2(P2/a2 2) [ l-2i/(k02r)-2/ (k02r)2]exp(-ika2r)
[0020] +A3[2i/(kpr)+2/(kpr)2]exp(-ikpr) (2a)
[0021] ?2=Αι(β2/αι2)[?/(k0ir)+l/(kair)2]exp(-ikair)+A2(P 2/a22)[i/(k02r) + l/ (k02r)2] exp(_ika2r)
[0022] +A3[l~i/(kpr)-l/(kpr)2]exp(-ikpr) (2b)
[0023] 式中系數(shù)Ai = (X2-X3) /(X2-X1); A2 = (X1-X3)/ (X1-X2); A3 = 1 (下同）
[0024] 應(yīng)力格林函數(shù)Tij表達式如下：
[0026] 其中：〇為單元應(yīng)力;η為單元法向量
[0027] 同理，另給出流體相對位移Wij、孔隙水壓力Pj格林函數(shù)表達式如下：
[0038] 2本單元格林函數(shù)解
[0039] 2.1集中力作用：
[0040] (1)固體骨架位移
[0041 ]由Gij(x，C) = 1^3^+(^2) γ i γ」]/4πμ；τ進行推導(dǎo)，其過程以水平方向力Fx作用 產(chǎn)生的X方向位移為例，為了方便，以離散單元中心為原點建立局部極坐標(biāo)系0-Θ，將水平力 F4V解為法向Fz，和切向Fx，，即Fz，= Fxnx
[0042] 在法向分解力Fz'作用下：
[0043] Gz,z,(xn,Ci) = [f2+(fi-f2) γζ· γz']/43iyr=f2/43iyr； γζ,=〇 (10)
[0044] 則法向分解力Fz，作用下產(chǎn)生的χ方向位移為：
[0046]在切向分解力Fx，作用下：
[0047] Gx'x' (χη, ξι) = [f2+(fi-f2) γ χ' γ χ' ]/4πμΓ= [f ι( l-cos20 )+f2( l+cos20) ]/8πμΓ ； γ χ· =-cos9
[0048] (12)
[0049] 貝IJ切向分解力Fx，作用下產(chǎn)生的χ方向位移為：
[0051]由公式（11)、（13)得出，水平方向力Fx作用在整體坐標(biāo)系下產(chǎn)生的χ方向的總位移 為：
[0053]則集中力作用(Fx，F(xiàn)y，F(xiàn)z)下的位移格林函數(shù)積分結(jié)果為：
[0055] 其中，Δ Si為離散單元面，δ。為Kronecker delta(當(dāng)i = j時，δ?」=1;當(dāng)i辛j時，δ。 =〇)，nj為尚散單兀法向量;式（15)中Fi、F2分別表不式⑵中fi、f2從〇SjR積分，R為尚散單兀 半徑，積分結(jié)果表達式如下：
[0063] (3)孔隙水壓力
[0064]由對稱性原理可知：
[0068] 當(dāng)^二纟:時，被積函數(shù)Tij具有奇異性其柯西主值為零，即當(dāng)離散單元為光滑圓面 時，牽引力格林函數(shù)張量Ti j ( Χη，ξ 1)為零，則：
[0069] tij(xn,Ci)=〇.55ij； (21)
[0070] 2.2流量源作用
[0071] (1)固體骨架位移：由式(5)、（7)可知，Pj和GiF之間的關(guān)系為:Pj = i wGiF由對稱性 可知，當(dāng)Χη=ξ?時：
[0073] (2)流體相對位移

[0085]所述散射波場通過施加在沉積盆地模型邊界上的虛擬波源與格林函數(shù)求得，具體 表達為：
[0090] 式中，Φ」，別為模型邊界固體骨架j方向虛擬波源和流量源;Gijjijjij和Pj 依次為集中力作用下的固體骨架位移函數(shù)、流體相對位移格林函數(shù)、牽引力格林函數(shù)和孔 壓格林函數(shù)，下標(biāo)j = l，2,3分別對應(yīng)虛擬波源在直角坐標(biāo)系下三個正交方向X，y，z;GlF， WiF，TiF和PF依次為流量源作用下的固體骨架位移函數(shù)、流體相對位移格林函數(shù)、牽引力格 林函數(shù)和孔壓格林函數(shù)。
[0091] (三)總波場計算
[0092] 根據(jù)場地邊界條件，即自由表面零應(yīng)力和耦合邊界位移、應(yīng)力連續(xù)條件構(gòu)造矩陣 方程組：[Η] [ Φ ] = [B]，并求解虛擬荷載密度[Φ ] = [H^B]，將所求得的虛擬波源和流量 源用于求解散射波場= 士⑷進而疊加自由波場得總波 場，即〉+ 5，其中自由波場表達式如下：
[0094] 式中，Φ和Φ為波勢函數(shù);kai、ka2和ke為波數(shù)；Θ為入射波和反射波在直角坐標(biāo)系下 與垂直方向的夾角；
[0095] 所述進行飽和格林函數(shù)的推導(dǎo)，采用間接邊界元法解決飽和局部場地的地震動模 擬，得出飽和局部場地的地震動模擬量化條件。
[0096] 本發(fā)明的效果是該方法提出了針對三維飽和場地對地震波散射的新的無奇異間 接邊界元法，并能很好的模擬固液動力耦合效應(yīng)、孔壓及流量的變化規(guī)律;為實際復(fù)雜工程 的模擬提供有效方法，為地震動響應(yīng)模擬提供理論依據(jù);能夠高效精確求解中小型復(fù)雜飽 和盆地地震動模擬和地震波散射模擬問題;為大規(guī)?；虼蟪叨热S飽和場地散射問題的快 速高效模擬奠定了基礎(chǔ)。該方法實現(xiàn)局部飽和不規(guī)則場地對地震波的散射模擬，以獲得該 地不同區(qū)域的地震動響應(yīng)，最終為城市地震區(qū)域化提供必要的理論依據(jù)。為飽和不均勻介 質(zhì)工程波動問題求解提供一種新的高效算法和實用程序，為沿海復(fù)雜飽和場地區(qū)域地震區(qū) 劃、結(jié)構(gòu)抗震設(shè)防提供理論依據(jù)，為實際復(fù)雜工程如地下隧道、復(fù)雜建筑結(jié)構(gòu)等的抗震模擬 帶來顯著的經(jīng)濟效益。
【附圖說明】
[0097] 圖1為本發(fā)明模擬方法實施流程圖；
[0098]圖2-a、圖2-b、圖2-c、圖2-d、圖2-e為本發(fā)明三維飽和盆地模型和邊界單元離散 圖；
[0099] 圖3-1為本發(fā)明離散單元整體坐標(biāo)系示意圖；
[0100] 圖3-2為本發(fā)明離散單元局部坐標(biāo)系示意圖；
[0101 ]圖4-1為本發(fā)明半球形沉積盆地三維計算模型圖；
[0102] 圖4-2為本發(fā)明半球形沉積盆地X〇Z截面；
[0103] 圖5為本發(fā)明平面P波垂直入射下沉積盆地地表位移幅值云圖；
[0104] 圖5-1為本發(fā)明孔隙率η = 0.1的X方向位移幅值；
[0105] 圖5-2為本發(fā)明孔隙率η = 0.3的X方向位移幅值；
[0106]圖5-3為本發(fā)明孔隙率η = 0.34的X方向位移幅值；
[0107] 圖5-4為本發(fā)明孔隙率η = 0.1的y方向位移幅值；
[0108] 圖5-5為本發(fā)明孔隙率η = 0.3的y方向位移幅值；
[0109 ]圖5-6為本發(fā)明孔隙率η = 0.34的y方向位移幅值；
[0110] 圖5-7為本發(fā)明孔隙率η = 0.1的z方向位移幅值；
[0111] 圖5-8為本發(fā)明孔隙率η = 0.3的ζ方向位移幅值；
[0112] 圖5-9為本發(fā)明孔隙率η = 0.34的ζ方向位移幅值。
【具體實施方式】
[0113] 結(jié)合附圖對本發(fā)明的三維飽和半空間中沉積場地地震動模擬方法加以說明。
[0114] 本發(fā)明的三維飽和半空間中沉積場地地震動模擬方法設(shè)計思想是基于結(jié)合邊界 元法，直接將邊界進行離散，使波源點與邊界點重合，根據(jù)邊界條件建立方程求解，將求得 的虛擬波源用于散射波場并疊加自由波場即得總波場，相比域離散、解析解等方法具有降 維、復(fù)雜模型適應(yīng)性好以及高精度等優(yōu)點。以飽和全空間中斜面圓盤荷載及流量源動力格 林函數(shù)為基本解，針對飽和兩相介質(zhì)發(fā)展一種新的無奇異間接邊界元法，并能很好的模擬 固液動力耦合效應(yīng)、孔壓及流量的變化規(guī)律。為實際復(fù)雜工程的模擬提供了一種有效方法， 為地震動響應(yīng)模擬提供理論依據(jù)。
[0115] 如圖1所示的模擬方法實施流程圖，本發(fā)明的三維飽和半空間中沉積場地地震動 模擬方法具體步驟為：
[0116] 1建立三維飽和場地模型：
[0117] 考慮平面波P或S和Rayleigh波的入射，基于單層位勢理論，在沉積盆地模型邊界 表面施加三個正交方向的虛擬荷載以構(gòu)造散射波場，采用邊界元方法實施中，需首先將邊 界表面離散為三角形或四邊形單元，為便于處理所述三角形或四邊形單元格林函數(shù)，采用 斜面圓盤均布荷載近似覆蓋在所述邊界表面的離散單元上，即分所述離散單元上虛擬荷載 對自身單元作用為本單元，所述離散單元上虛擬荷載對自身單元之外的離散單元作用為非 本單元，推導(dǎo)出飽和介質(zhì)中三個方向的力F x，F(xiàn)y，F(xiàn)z集中力和流量元Q作用下的格林函數(shù)，進 而根據(jù)所述邊界表面的透水和不透水邊界條件構(gòu)造矩陣方程組：
[0118] [Η] [ Φ ] = [B]，求解得到虛擬荷載密度：[Φ ] = [Η]-KB]
[0119]其中，[Η]為格林函數(shù)影響矩陣，[B]為邊界條件解。
[0120] 如圖2-a、圖2-b、圖2-c、圖2-d、圖2-e所示的三維飽和盆地模型和邊界單元離散 圖。根據(jù)透水和不透水邊界條件建立方程組，求解得到虛擬荷載密度。散射波場由各個單元 上的虛擬荷載的綜合作用而得，最后將散射波場和自由波場疊加即得到總波場。由于虛擬 荷載直接施加在邊界面上，本方法對邊界形狀具有更好的適應(yīng)性，可以解決任意形狀的局 部場地。
[0121] 1 Biot飽和兩相介質(zhì)理論
[0122] 基于Biot理論，設(shè)土體骨架位移和流體相對于土體骨架的位移分別為udPWl(i = X，y，Z)，則均勻飽和介質(zhì)本構(gòu)關(guān)系表達式：
[0123] 〇ij = 2yeij+A5ije-a5ijp； i , j = x,z (la)
[0124] p = -aMe-Mwi,i (lb)
[0125] 與u，w相關(guān)的多孔介質(zhì)運動方程表達式：
[0126] + (A + α~Μ + μ)ιι]β + aMwjfi = put + (2a)
[0127] α/V/i^ ；/ + Mw. μ - pfii. +mn\ + b\\\ (2b)
[0128] 其中，Pi (快波）、P2 (慢波）的勢函數(shù)Φ l、Φ 2和S波勢函數(shù)Φ、x分別滿足以下波動方 程：
[0131]基于Helmholtz矢量分解原理，在柱坐標(biāo)系下土體骨架位移u和流體相對骨架位移 w為：
[0136] 其中，Φ i、Φ 2、Ψ和X為孔隙流體相應(yīng)勢函數(shù)，滿足下列等式：
[0137] Φ1 = χ1 Φ 1; Φ2 = χ2 φ 2; ψ =χ3φ; X =χ3χ (5)
[0138] 參數(shù) XiAA 為：
[0146] 另外定義飽和土中的剪切波波數(shù)為:? =(P>3 + 1)歲。
[0147] 2三維飽和格林函數(shù)推導(dǎo)
[0148] 2.1非本單元格林函數(shù)解
[0149] 2.1.1集中力作用
[0150] 類似于三維彈性半空間問題，在均勻各向同性彈性固體介質(zhì)中，位移格林函數(shù)可 表達為：
[0151] Gij(x,C) = [f25ij+(fi-f2) γ i γ j]/43iyr (11)
[0152] 其中，yj = (XjU/r，r表示波源點（ξχ,ξιξΟ與接收點（X1，X2， X3)之間的距離，即 1'2=(叉1-|1)2+(叉2-|2) 2+(叉3-|3)2,31」表示1(1'〇116。1^1(16]^，4表示拉梅常數(shù)；1^=〇/0，1^1 =?/ <11，1^2=(〇/(12分別表示5，？ 1(快波），？2(慢波)波數(shù);{3，(11，(12表示其各自對應(yīng)波速 (^1和 f2定義為：
[0153] ?>ι = Αι(β2/α12) [l-2i/(kair)-2/(kair)2]exp(-ikair)+A2(02/a2 2)[卜2i/(ka2r)-2/ (ka2r)2]exp(-ika2r)
[0154] +As[2i/(k^r)+2/(k^r)2]exp(-ik^r) (12a)
[0155] f2=Ai(02/a12)[i/(kair)+l/(kair) 2]exp(-ikair)+A2(02/a22)[i/(ka2r) + l/(ka2r)2] exp(-ika2r)
[0156] +As[l-i/(k^r)-l/(k^r)2]exp(-ik^r) (12b)
[0157] 式中系數(shù)六1=(乂2-乂3)/(乂2-乂1);六2=(乂1-乂3)/(>1-乂2);六3=1(下同）
[0158] 應(yīng)力格林函數(shù)Tij由本構(gòu)關(guān)系給出：
[0161]同理，另給出其他格林函數(shù)表達式，即流體相對位移Wij、孔隙水壓力Pj表達式如 下：
[0165] f ,=χιΑι(β2/αι2) [ l_2i/(kair)-2/(kair)2]exp(-ikair )+Χ2Α2(β2/α22) [ l-2i/(ka 2r)-2/ (ka2r)2]exp(-ika2r)
[0166] +X3As[2i/(k^r)+2/(k^r)2]exp(-ik^r) (16a)
[0167] f2w=xiAi(P2/ai2) [i/(kair)+l/ (kair)2]exp(-ikair)
[0168] +Χ2Α2(β2/α22) [i/(ka2r)+l/(ka2r)2]exp(-ika2r)
[0169] +X3A3[l-i/(ksvr)-l/(ksvr) 2]exp(-iksvr) (16b)
[0170] 2.1.2流量源作用
[0171 ] GiF(x,C)=-iB γ i[exp(-ikair)(ikair+l)/r-exp(-ika2r)(ika2r+l)/r]/4Ji〇r (17)
[0172] ffiF(x,C)=_iBy i[xiexp(-ikair)(ikair+l)/r-X2exp(-ika2r)(ik02r+l)/r]/4Ji〇r (18)
[0173] (λ%^) = iBM[k^(a -f /^^xpi-ih^r) / r - kl2(a -f χ2 )A2 Qxp(-ikalr) / r] / Απω ( 19)
[0174] 其中：B = 1/(Xi_X2)，系數(shù)Xi，X2，X3見式(6)。
[0175] 2.2本單元格林函數(shù)解
[0176] 2.2.1集中力作用：
[0177] (1)固體骨架位移
[0178]如圖3-1、圖3-2所示的離散單元局部極坐標(biāo)系示意圖，由Gij(x，ξ) = [f2δij+(f1-f2)γlγj]/4πμr進行推導(dǎo)，其過程以水平方向力F x作用產(chǎn)生的X方向位移為例，為了方便， 以離散單元中心為原點建立局部極坐標(biāo)系〇-θ，將水平力Fx分解為法向Fz，和切向Fx，，gPF z， = Fxnx，
[0179]在法向分解力Fz，作用下：
[0180] Gz,z,(xn,Ci) = [f2+(fi-f2) J z γ z']/4πμΓ = f2/4jiyr ； γζ'=0 (20)
[0181] 則法向分解力Fz，作用下產(chǎn)生的x方向位移為：
[0183]在切向分解力Fx，作用下：
[0184] Gx'x' (χη, ξι) = [f2+(fi-f2) γ x' γ x' ]/4πμΓ= [f i( l-cos20 )+f2( l+cos20) ]/8πμΓ ； γ x' =-cos9
[0185] (22)
[0186] 則切向分解力Fx，作用下產(chǎn)生的x方向位移為：
[0188] 由此可得，水平方向力Fx作用在整體坐標(biāo)系下產(chǎn)生的X方向的總位移為：
[0189] = jj' [(./； + Id HA -Ιι)ηΛ(24)
[0190] 則集中力三個方向Fx，F(xiàn)y，F(xiàn)z作用下的位移格林函數(shù)積分結(jié)果為：（其他方向作用下 的推導(dǎo)與水平方向力作用類似）
[0192]其中，Δ Si為離散單元面，δ。為Kronecker delta(當(dāng)i = j時，5ij = 1;當(dāng)j時，δ。 =〇)，nj為禺散單兀法向量。式(25)中Fi、F2分別表亦式（12)中fi、f2從0到R積分，R為禺散單 元半徑。積分結(jié)果表達式如下：

[0205] 當(dāng)知二纟:時，被積函數(shù)Tij具有奇異性其柯西主值為零，即當(dāng)離散單元為光滑圓面 時，牽引力格林函數(shù)張量Ti j ( Χη，ξ 1)為零，則：
[0206] tij(xn,Ci)=0.55ij (31)
[0207] 2.2.2流量源作用
[0208] (1)固體骨架位移：由式（15)、（17)可知，Pj和GiF之間的關(guān)系為:Pj = ic〇GiF由對稱 性可知，當(dāng)Χη=ξ?時：
[0223]至此，本發(fā)明中散射波場即可通過施加在邊界上的虛擬波源與格林函數(shù)求得，其 具體表達為：

[0228]式中，Φ」，別為模型邊界固體骨架j方向虛擬波源和流量源;Gijijjij和Pj 依次為集中力作用下的固體骨架位移函數(shù)、流體相對位移格林函數(shù)、牽引力格林函數(shù)和孔 壓格林函數(shù)，下標(biāo)j = 1，2，3分別對應(yīng)虛擬波源三個正交方向X，y，z ;GiF，WiF，TiF和PF依次為 流量源作用下的固體骨架位移函數(shù)、流體相對位移格林函數(shù)、牽引力格林函數(shù)和孔壓格林 函數(shù)。
[0229]最后，根據(jù)場地邊界條件，即自由表面零應(yīng)力和耦合邊界位移、應(yīng)力連續(xù)條件構(gòu)造 矩陣方程組：[Η] [ Φ ] = [B]，并求解虛擬荷載密度[Φ ] = [ΗΓΚΒ]，將所求得的虛擬波源和 流量源用于求解散射波場，疊加自由波場即可得總波場。
[0230] 三維飽和場地地震響應(yīng)分析，參考沿海地區(qū)實際沉積盆地復(fù)雜的地質(zhì)環(huán)境、幾何 和材料特性，針對復(fù)雜飽和沉積盆地對地震波的三維散射，基于本發(fā)明三維飽和ΙΒΕΜ，進行 地震波散射模擬。沉積盆地地形作為常見的局部場地之一，研究在地震作用下的場地反應(yīng) 具有十分重要的理論和現(xiàn)實意義，可以為建筑結(jié)構(gòu)或橋隧工程的抗震設(shè)計中提供科學(xué)定量 的參數(shù)分析，提高實際工程的抗震性能。
[0231] 如圖4-1、4_2所示，以飽和半空間中半球形沉積盆地為例，采用本發(fā)明三維飽和 ΙΒΕΜ求解了其對縱波（Ρ波）的散射。半空間飽和介質(zhì)：臨近孔隙率= 0.36，泊松比υ = 0.25，材料阻尼ζ = 0.002，土體骨架體積模量Kgs = 36000Mpa，土體顆粒密度pgs = 2650kg/m3， 土體臨界體積模量Kcr = 200Mpa，流體密度pf = 1000kg/m3，流體體積模量Kf = 2000Mpa;沉積 內(nèi)飽和介質(zhì)：臨界孔隙率nCT = 0.36，泊松比υ = 〇. 25，材料阻尼ζ = 0.002，土體骨架體積模量 Kgv = 9000Mpa，土體顆粒密度pgv = 2650kg/m3，土體臨界體積模量Kcr = 50Mpa，流體密度Pf = 1000kg/m3，流體體積模量Kf = 2000Mpa。入射波頻率η = 0 · 5，孔隙率η = 0 · 1，0 · 3和0 · 34。
[0232] 如圖5-1-圖5-9所示的沉積盆地地表位移幅值云圖，通過大量數(shù)值模擬結(jié)果表 明：本發(fā)明的模擬方法能夠高效精確求解三維飽和不規(guī)則場地對地震波的散射問題，相比 單相的彈性半空間，飽和半空間場地對地震波的散射特性存在顯著差異，且對于不同入射 頻率的地震波其場地地震響應(yīng)變化劇烈，對于無限域問題，能夠很好的避免由于場地截斷 引起的誤差。同時，隨孔隙率的增加，盆地地表位移幅值放大效應(yīng)增強。
[0233] 上面對本發(fā)明的較佳實施方式作了詳細說明，但是本發(fā)明并不限于上述實施方 式，在本領(lǐng)域的普通技術(shù)人員所具備的知識范圍內(nèi)，還可以在不脫離本發(fā)明宗旨的前提下 做出各種變化均為本權(quán)利要求保護范圍。
【主權(quán)項】
1. 一種三維飽和半空間中沉積場地地震動模擬方法，該方法是在建立三維飽和沉積盆 地模型基礎(chǔ)上，通過由飽和土中的土體骨架和流體之間的相對關(guān)系，推導(dǎo)飽和格林函數(shù)，采 用間接邊界元法解決飽和局部場地的地震動模擬，該方法包括有以下步驟:量化 (一) 建立三維飽和場地模型： 考慮平面波P或S和Rayleigh波的入射，基于單層位勢理論，在沉積盆地模型邊界表面 施加三個正交方向的虛擬荷載以構(gòu)造散射波場，采用邊界元方法實施中，需首先將邊界表 面離散為三角形或四邊形單元，為便于處理所述三角形或四邊形單元格林函數(shù)，采用斜面 圓盤均布荷載近似覆蓋在所述邊界表面的離散單元上，即分所述離散單元上虛擬荷載對自 身單元作用為本單元，所述離散單元上虛擬荷載對自身單元之外的離散單元作用為非本單 元，推導(dǎo)出飽和介質(zhì)中三個方向的力F x，F(xiàn)y，F(xiàn)z集中力和流量元(Q)作用下的格林函數(shù)，進而 根據(jù)所述邊界表面的透水和不透水邊界條件構(gòu)造矩陣方程組： [Η] [ Φ ] = [B]，求解得到虛擬荷載密度：[Φ ] = [Η]ΛΒ] 其中，[Η]為格林函數(shù)影響矩陣，[Β]為邊界條件解 最后將邊界表面上各個單元上的虛擬荷載的綜合作用而得散射波場和不含局部不均 勻場地時彈性波入射下的半空間自由波場疊加即得到總波場；由于虛擬荷載直接施加在邊 界表面上，能夠解決多種形狀的飽和局部場地； (二) 三維飽和格林函數(shù)推導(dǎo)： 基于Biot飽和兩相介質(zhì)理論，推導(dǎo)飽和介質(zhì)中三個方向的力Fx，F(xiàn)y，F(xiàn)z集中力和流量元 (Q)作用下的格林函數(shù)，格林函數(shù)解如下： 1非本單元格林函數(shù)解 1.1集中力作用 類似于三維彈性半空間問題，在均勻各向同性彈性固體介質(zhì)中，位移格林函數(shù)表達為： Gij(x,C) = [f25ij+(fi-f2) γ i γ j]/43iyr (1) 其中，丫 j = (xjU/r，r表示波源點（ξι,ξ2,ξ3)與接收點(xi,X2,X3)之間的距離，即r2 = (χι-ξι)2+(Χ2-ξ2)2+(Χ3_ξ3)2，δ?』表示Kronecker delta，y表示拉梅常數(shù)；ksv= ω/β，1?α?= ω/ <11，1^2=(〇/(12分別表示5，？1(快波），？ 2(慢波)波數(shù);{3，(11，(12表示其各自對應(yīng)波速，&和5定 義為： ?'ι = Αι(β2/α12) [ l-2i/(kair)-2/(kair)2]exp(-ikai;r)+A2(i52/a2 2) [l-2i/(ka2r)-2/(ka2r )2]exp(-ik〇2r) +A3[2i/(kpr)+2/(kpr)2]exp(-ikpr) (2a) ?'2 = Αι(β2/α12) [i/(kair)+l/(kair)2]exp(-ikai;r)+A2(02/a2 2) [i/(ka2r) + l/(ka2r)2]exp (_ika2r) +A3[ l~i/(kpr)-l/ (kpr)2]exp(-ikpr) (2b) 式中系數(shù)Ai = (X2-X3)/ (X2-X1); A2 = (X1-X3) /(X1-X2); A3 = 1 (下同） 應(yīng)力格林函數(shù)Tij表達式如下： (3) (4) 其中：σ為單元應(yīng)力;η為單元法向量 同理，另給出流體相對位移Wij、孔隙水壓力Pj格林函數(shù)表達式如下： Pj(χ,ξ) = MYj[k1al{a + ^)4 Qxp(-ikair)(ikalr +1)/ r + kl2(a + χ2)A2 exp(-ika2r)(ika2r +1)/r]/ Απμι-kf (5) 其中Γ和/7為：1.2流量源作用 Gif(x,ξ) = -iB γ i[exp(-ik〇ir) (ikair+l)/r-exp(-ika2r) (ik〇2r+l)/r]/4jr ωr (7) WiF(x,ξ) = _iB γ i[xiexp(_ikair) (ikair+1 )/r_X2exp(_ika2r) (ika2r+l)/r]/4iT ωr (8)其中：B=1/(xi-x2),系數(shù)i 2本單元格林函數(shù)解2.1集中力作用： (1)固體骨架位移 由61辦，|) = [0"+(^2)丫1丫」]/43^進行推導(dǎo)，其過程以水平方向力?\作用產(chǎn)生的 X方向位移為例，為了方便，以離散單元中心為原點建立局部極坐標(biāo)系〇-θ，將水平力Fx分解 為法向Fz，和切向Fx，，即Fz，=Fxn x，F(xiàn),. =： F^\-n；； 在法向分解力Fz，作用下： Gz,z,(xn,Ci) = [f2+(fi-f2) γζ· γ z']/43iyr = f2/43Tyr； γζ,=0 (10) 則法向分解力F7.，作用下產(chǎn)牛的χ方向位移為：在切向分解力Fx，作用下： Gx'x'(χη,ξι) = [f2+(fi-f2) γχ· yx,]/43Tyr=[fi(l-cos20)+f2(l+cos20)]/83ryr； γχ,=- cos9 (12) 則切向分解力Fx，作用下產(chǎn)生的χ方向位移為：由公式(11 )、（ 13)得出，水平方向力Fx作用在整體坐標(biāo)系下產(chǎn)生的X方向的總位移為：則集中力作用（Fx,Fv,F,)下的份務(wù)格林函數(shù)積分結(jié)果為：其中，Δ Si為離散單元面，δ。為Kronecker delta(當(dāng)i = j時，δ?」= 1;當(dāng)i辛j時，δ?」= 〇)， nj為咼散單兀法向量;式（15)中Fi、F2分別表不式⑵中fi、f2從〇SjR積分，R為尚散單兀半徑， 積分結(jié)果表達式如下： 、 ·>.·[/ ，-a 、 J V， V t冰丄丨'14; L·. V l> f 7 \ - f (2) 流體相對位移其中：(3) 孔隙水壓力 由對稱性原理可知：(4) 牽引力當(dāng)χη=ξι時，被積函數(shù)Tu具有奇異性其柯西主值為零，即當(dāng)離散單元為光滑圓面時，牽 弓丨力格林函數(shù)張量Ti j ( Χη，ξ 1)為零，則： tij(xn,Ci)=0.55ij； (21) 2.2流量源作用 (1) 固體骨架位移：由式(5)、（7)可知七和61「之間的關(guān)系為七=1〇61「由對稱性可知， 當(dāng)Χη=ξ?時：1 流體相對位移其中，0' :-:(m如丨μ ^χη = ξι時，WiF(xn，li)為零，貝丨J:所述散射波場通過施加在沉積盆地模型邊界上的虛擬波源與格林函數(shù)求得，具體表達 為：式中，Φ」，(^分別為模型邊界固體骨架j方向虛擬波源和流量源;Gijjijjij和Pj依次 為集中力作用下的固體骨架位移函數(shù)、流體相對位移格林函數(shù)、牽引力格林函數(shù)和孔壓格 林函數(shù)，下標(biāo)j = 1，2，3分別對應(yīng)虛擬波源在直角坐標(biāo)系下三個正交方向X，y，z ; GiF，WiF，TiF 和Pf依次為流量源作用下的固體骨架位移函數(shù)、流體相對位移格林函數(shù)、牽引力格林函數(shù) 和孔壓格林函數(shù)； (三)總波場計算 根據(jù)場地邊界條件，即自由表面零應(yīng)力和耦合邊界位移、應(yīng)力連續(xù)條件構(gòu)造矩陣方程 組：[Η] [ Φ ] = [B]，并求解虛擬荷載密度[Φ ] = [H]-乜]，將所求得的虛擬波源和流量源用 于求解散射波場+，進而疊加自由波場得總波場，即 (32) 〇(，) +?f):，其中自由波場表達式如下： 式中，Φ和Φ為波勢函數(shù);kal、ka2和kf!為波數(shù)；Θ為入射波和反射波在直角坐標(biāo)系下與垂 直方向的夾角； 所述進行飽和格林函數(shù)的推導(dǎo)，采用間接邊界元法解決飽和局部場地的地震動模擬， 得出飽和局部場地的地震動模擬量化條件。
【文檔編號】G06T17/05GK105869208SQ201610184315
【公開日】2016年8月17日
【申請日】2016年3月28日
【發(fā)明人】劉中憲, 劉蕾, 王冬, 王海良, 尚策
【申請人】天津城建大學(xué), 中煤科工集團北京華宇工程有限公司平頂山分公司