流固耦合地震波動問題的顯式譜元模擬方法1)

孔曦駿 邢浩潔 李鴻晶,2)

* (南京工業大學工程力學研究所,南京 211816)

? (中國地震局地球物理研究所,北京 100081)

引言

流固耦合地震波動是地震工程領域中近年來備受重視的一個問題.處理流固耦合問題的一般性做法是認為流體和固體兩種介質的耦合作用僅發生在交界面上,因而只需引入界面協調條件以實現波動能量的交換.這種傳統分析方法要求流體和固體中的波動以不同的方程進行描述,建立的數值格式不統一,界面耦合處理亦很復雜,應用上十分不方便.陳少林等[1-2]利用兩相介質理論發展了流固耦合問題的統一計算框架,將流體和固體分別視為飽和多孔介質退化后的特殊形態,從而通過統一的數值格式實現對流固耦合波動的模擬.流固耦合問題統一計算框架避免了不同模塊之間的交互操作,易于實現并行計算,為構建高效率的時空解耦顯式算法以實現大規模的流固耦合波動模擬提供了可能.

流固耦合地震波動模擬方法主要有解析法、半解析法、邊界元法、時域及頻域有限元法、有限差分法等.杜修力等[3]、趙成剛等[4]和王進廷等[5-6]基于集中質量有限元和外推人工邊界發展出一套適用于重力壩/拱壩-庫水-飽和場地-基巖系統地震反應分析的數值模擬方法,全面探討了流固耦合界面處理、顯式時間積分格式、地震動輸入、大型復雜系統建模等問題.李偉華[7]和陳少林等[8-9]對飽和兩相介質場地中的波傳播與散射理論、顯式有限元方法以及人工邊界條件進行了進一步研究.李亮等[10-11]、宋佳等[12-13]基于u-p格式的飽和波動方程研究了時域全顯式有限元模擬方法以及飽和波動的應力型人工邊界條件.梁建文等[14]、巴振寧等[15]、劉中憲等[16]基于多極子快速間接邊界元法研究了飽和層狀場地對地震波的二維和三維散射.Wang 等[17-19]、Zhao等[20-21]采用半解析法、頻域有限元法研究了海水-飽和海床-基巖系統的地震反應以及單樁海上風機、防波堤結構的動力反應.Chen 等[22-23]使用頻域有限元法研究了考慮土-結構-海水相互作用時單樁基礎以及海堤的頻域動力響應.Sun 等[24]通過曲線網格有限差分法研究了考慮海水作用時三維海底場地的地震反應特征.Liu 等[25]、Bao 等[26]、寶鑫等[27-28]基于可壓縮流體-彈性固體(考慮非線性)的流固耦合模式,將自主開發的固體介質黏彈性人工邊界、流體介質動力人工邊界、邊界子結構/內部子結構的地震動輸入方法、成層介質側邊界地震動輸入的一維化時域方法等技術與商業有限元軟件相結合,實現了島礁-海水系統地震反應的有限元分析.陳國興等[29-30]對典型的實際海域島礁場地以及海峽場地的二維剖面進行了非線性地震反應分析,探討了海域場地條件對地震動參數的影響規律.

現有方法中,解析法/半解析法通常只適用于比較規則或由特定函數描述的幾何形狀,同時由于計算過程中高階級數項數值很小容易導致溢出計算機精度范圍,使得這類方法能夠模擬的頻率上限不高.有限差分法目前常用高階差分格式來提高模擬波場的精度并且能夠模擬不均勻介質,但是其處理復雜地表地形及介質交界面時往往精度不高并且容易引起穩定性問題,另外它在模擬面波以及界面波方面的效果較差.邊界元法能夠模擬很寬的頻帶范圍,且近年來發展出多極子技術來提高計算速度,已被成功地應用于大型河谷飽和場地以及城市尺度沉積盆地的地震動模擬,但是其計算模型主要為單一介質的不規則區域或橫觀各向同性的水平成層及圓弧成層場地.有限元法最大優勢在于通過靈活的單元劃分來實現對各種復雜幾何形狀的建模,以及從建立單元到集成為整體計算模型過程的規整性,為開發大型計算程序帶來了便利,在應用于場地地震波動模擬時又開發出了集中質量有限元技術以提高計算效率,不過,該方法目前主要使用基于一階插值函數的四邊形單元(二維)或六面體單元(三維),這使得它模擬波場的精度比較有限,高頻段模擬結果往往并不準確.

譜元法[31-34]可以看作一種特殊的高階有限元法,其優點在于使用高階單元(如四階以上的單元)進行波動模擬并且從理論上嚴格地導出集中質量矩陣,成功地解決了目前基于低階單元(常用一階單元)的集中質量有限元法模擬波場的精度不高以及現有基于經驗的集中質量方法缺乏嚴格數學基礎的問題.Duczek 等[35]證明了譜元法基于節點積分直接導出的集中質量矩陣,與傳統基于行和集中、對角元素放大等經驗方法得出的集中質量矩陣結果一致.近二三十年來的研究和應用已較為清楚地表明,譜元法是有限元法用于各類波動問題數值模擬時一種最為高效(既能高精度地模擬波場又能靈活地對復雜幾何形狀進行建模且大規模問題的計算量可接受)的形式.

本文以飽和兩相介質波動的Biot 方程及其向理想流體(孔隙率為1 情形)、彈性固體(孔隙率為0 情形)的兩種退化形式方程[1-2,32,36]為基礎,采用高精度譜元法與多次透射公式人工邊界條件[37-40]相結合,發展流固耦合地震波動問題的高效數值模擬方法.通過與傳遞函數法[41]以及集中質量有限元法的結果進行比較,驗證本文方法的高效性和地震波動模擬結果的寬頻帶特性.

1 流固耦合問題的控制方程

考慮如圖1 所示流固耦合問題,該問題模型由理想流體、飽和多孔介質和彈性固體三種介質組成,亦可以退化為只含有兩種介質之間耦合的情形.

圖1 流固耦合地震波動問題示意圖Fig.1 Schematic diagram of seismic wave motion problem in a fluidporoelastic-solid model

本文流固耦合地震波動問題的控制方程由Biot 所提出的飽和多孔介質波動方程及其向理想流體、彈性固體的兩種退化形式共同組成.這樣組合的優點是飽和波動的求解模式可以直接應用于理想流體和彈性固體,并且飽和多孔介質之間的界面連續條件亦可以轉化為任意兩種介質界面情形,接下來分別進行介紹.

考慮二維情形,當場變量為固相骨架和液相流體的位移時,Biot 流體飽和多孔介質波動方程為

式中,u,U分別為固相骨架位移、液相流體位移的列向量,u=(ux,uy)T,U=(Ux,Uy)T,變量上方一個、二個圓點分別表示對時間求一階、二階導數,上標“T”表示矩陣或向量轉置;Ls和Lw為微分矩陣,D1,D2和D3為介質的力學參數,Lw=[?/?x,?/?y],D2=Q,D3=R

這里β為孔隙率,ρs為固相骨架密度,ρw為液相流體密度,參數b=β2μ0/k0,μ0為孔隙流體的動黏度系數,k0為達西滲透系數;A,G,Q,R為Biot 常數,可通過如下表達式計算

其中λ和μ為固相骨架的拉梅常數,和M為表征飽和多孔介質壓縮性的常數.在實際應用中,相關常數可以通過試驗來測出.

Biot 方程(1)當孔隙率為1 時退化為理想流體方程.此時孔隙率β=1,以及常數α=1,且由于不存在固相骨架,動黏度系數為零,導致參數b=0,固相模量亦置為0.于是Biot 常數為:A=0,G=0,Q=0,R=M.方程(1)改寫為α

Biot 方程(1)當孔隙率為0 時退化為彈性固體的波動方程.此時孔隙率β=0,以及常數α=0,且由于不存在固相骨架,動黏度系數為零,導致參數b=0.于是Biot 常數為A=λ,G=μ,Q=0,R=0.方程(1)改寫為

類似上述由Biot 方程導出理想流體和彈性固體方程的做法,不同介質交界面連續條件同樣可以由兩種不同力學參數飽和多孔介質的交界面這種一般情形,分別退化到理想流體–飽和多孔介質、飽和多孔介質-彈性固體、理想流體–彈性固體三種交界面情形,如圖2 所示.

圖2 流固耦合問題的不同介質交界面Fig.2 Different interfaces in the fluid-poroelastic-solid problem

兩種不同力學參數飽和多孔介質交界面的連續條件為:

(a) 交界面兩側介質的法向總應力相等;

(b) 交界面兩側固相骨架的切向應力相等;

(c) 交界面兩側液相流體的壓強相等;

(d) 交界面兩側固相骨架的法向、切向位移相等;

(e) 交界面兩側的法向運動連續.

上述界面連續條件可用表達式分別表示為

式中,上標“A”和“B”指示交界面兩側的不同介質,下標“N”和“T”分別指示變量的法向分量和切向分量;σ為固相骨架應力,τ為平均孔壓,τ=-βP,β為孔隙率,P為孔隙流體壓力.

式(6)所給出的飽和多孔介質A-飽和多孔介質B 這種一般情形下的界面連續條件向理想流體-飽和多孔介質、飽和多孔介質-彈性固體、理想流體-彈性固體三種具體情形下交界面連續條件的退化,可按照由方程(1)退化到方程(4)或方程(5)的基本規律以及由圖2 所給出的參數取值方式進行.不過,在退化過程中式(6)的一些項需要隨著固相或液相成分的消失而舍去,相關細節可參見文獻[1-2,5,36],本文不再贅述.

綜上,內域波動方程式(1)、式(4b)和(5a)和界面連續條件方程(6)組成了本文流固耦合地震波動問題的控制方程.

2 流固耦合地震波動的譜元分析方法

為提高流固耦合地震波動問題的求解精度和計算效率,不同于現有的集中質量有限元法,本文采用能夠天然地導出集中質量高階單元的譜元法來求解方程式(1)、式(4b)、式(5a)和式(6)組成的系統在地震作用下的動力反應.譜元法可以看作一種特殊的高階有限單元,因此其基本步驟與有限元法一致,區別僅在于它采用的是基于譜插值方式構建的高階單元.對上述流固耦合方程系統的空間離散,由于退化關系的存在,使得只需詳細探討對于飽和兩項介質波動Biot 方程(1)和一般情形界面連續條件方程(6)的譜元離散,然后根據退化規則即可退化到其他情形.本節首先論述內域波動方程的譜元離散及其時域逐步積分格式,然后在此基礎上介紹界面連續條件的離散化表示,最后探討人工邊界條件與地震動輸入,這樣形成一套流固耦合地震波動的譜元分析方法.

2.1 譜元離散方程

首先考慮飽和波動方程(1)的譜元離散.與有限元法類似,譜元離散亦需在原方程的等效積分弱形式上進行.對方程(1)的每一項乘以任意的測試函數v或V,對空間變量在全域上進行積分;對于方程右邊關于空間導數的各項進行分部積分,利用格林公式將域積分轉化為邊界積分,引入自由地表條件使得邊界積分項等于零而從方程中消失.得到等效積分弱形式方程如下

在弱形式方程(7)中,將計算模型Ω劃分為各個譜單元Ωe的組合,則全域積分轉換為各個譜單元積分的和.將各個譜單元的積分、微分運算轉化到[-1,1]×[-1,1]的標準參考單元Λ上進行.

根據Galerkin 離散化法則,對以單元形式表示的弱形式方程進行數值離散,單元內的各個場變量和測試函數均采用統一的譜插值模式來構建單元插值函數.譜單元的插值函數在標準參考坐標Λ(ξ,η)∈[-1,1]×[-1,1]內進行定義,其中單元場變量的插值函數表達式為

式中,ux,i,uy,i,Ux,i,Uy,i(i=1,2,···,ng)為待求解的離散網格節點上的場函數值,Ni(ξ,η) 為單元內第i個節點的形函數,ng為譜單元的節點數目.

形函數表達式為

式中,ngllξ和ngllη分別表示ξ,η方向的單元節點數目.lm(ξ),ln(η)均為Lagrange 插值基函數,ξj,ηj為Gauss-Lobatto-Legendre (簡稱GLL)高精度數值積分節點[35,39].本文采用常用的四階譜單元,其GLL 節點為[-1,-0.654 653 670 7,0,0.654 653 670 7,1].

根據上述譜單元場函數的定義,可以完成以單元形式表示的弱形式方程的數值離散,然后將各個單元離散方程集成為總體離散方程.于是,得到飽和兩項介質波動Biot 方程(1)采用譜元法進行數值離散的運動方程為

式中,u和U為單元節點位移的列向量,Ms為固相質量矩陣,Mw為液相質量矩陣,C為流固耦合阻尼矩陣(用于計算固液兩相之間的黏性阻力),Kss為固相剛度矩陣,Kww為液相剛度矩陣,Ksw和Kws為代表固液兩相之間本構力的耦合矩陣;Fs為作用在固相上的外荷載矩陣,Fw為作用在液相上的外荷載矩陣,FRs和FRw分別為單元之間作用在固相和液相上的相互作用力,主要體現在不同介質的分界面上,若單元在同一介質內時,該力可被忽略.上述各個整體矩陣由相應的單元矩陣集成得到

式中,N為形函數矩陣,J=J(x,y;ξ,η)為雅可比矩陣,|J|為雅可比行列式,矩陣Bs=LsN,矩陣Bw=LwN,為由方向導數組成的矩陣,n為沿外邊界外法線的方向矢量,相關表達式為

根據第1 節所揭示的波動微分方程退化關系,按照對應的參數取值規則,可以方便地由這里的飽和波動譜元離散運動方程(10)分別退化得到理想流體、彈性固體的離散運動方程.因此在計算程序中,只需適當地調整飽和多孔介質的單元特性矩陣及總體特性矩陣中的相關參數取值,就可以得到理想流體、彈性固體的單元特性矩陣及總體特性矩陣.這種統一的處理方式首先由陳少林等[1-2,36]提出,使得流固耦合問題的離散化建模過程大為簡化.

2.2 時域逐步積分

傳統有限元法中的總體質量矩陣Ms,Mw和流固耦合阻尼矩陣C均為非對角矩陣,顯式的時步積分格式需要對質量矩陣(Ms,Mw)進行人為的行和集中等措施才能實現,但流固耦合阻尼矩陣C通常不做處理,依然為非對角矩陣.因此,基于集中質量有限元的流固耦合地震波動問題一般選取具有一階計算精度的中心差分-單邊差分法或者形式繁瑣但具有二階計算精度的中心差分-Newmark 常平均加速度法.

譜元法基于GLL 積分直接導出集中質量矩陣.文獻[35]從理論角度嚴格證明了譜元質量模型與傳統基于行和集中、對角元素放大等經驗方法得出的集中質量矩陣在本質上是一致的,但具有堅實的數學基礎.然而很少有研究關注到Biot 方程中的流固耦合阻尼矩陣C在應用譜元法進行GLL 積分后同樣為對角矩陣,這一點可以從式11(a)、式11(b)與式11(c)的對比中直接看出: 質量矩陣和阻尼矩陣的積分公式只有系數上的區別,進行GLL 數值積分后皆為對角矩陣.

如何理解流固耦合阻尼矩陣C為對角矩陣的物理本質對波場問題的數值計算具有重要意義.有限元的行和集中質量矩陣被認為是忽略了單元內加速度的變化,即假定單元內慣性力為常量.基于GLL 積分的流固耦合阻尼矩陣C為對角矩陣,則可以類似地理解為忽略了單元內固、液兩相速度之差的變化,即假定單元內由兩相速度差導致的黏性阻力為常量.

考慮到本文中總體質量矩陣Ms,Mw和流固耦合阻尼矩陣C均為對角矩陣,可以采用計算精度為二階且形式簡潔的中心差分法來處理譜元離散方程(10)中的時間導數

其中,上標p表示t=pΔt時刻.將式(13)代入譜元離散方程(10),得到時域逐步積分格式

其中,u和U表示被計算節點的固、液相位移;ue和Ue表示與該計算節點相關單元內節點的固、液相位移.時域逐步積分格式(14)是關于整個譜元離散模型的時間步進格式.式中含有模型的整體特性矩陣,實際編程時這樣做很不經濟,且對于大規模問題往往難以實現.考慮到這里總體質量矩陣Ms,Mw和流固耦合阻尼矩陣C均為對角矩陣,所以式(14)可以寫成關于各個節點運動ux,uy,Ux,Uy的完全顯式的時間積分格式.以編號為k的節點固相x向位移ux為例,其顯式逐步積分格式可表示為

其中,nc為與節點k相關的節點的數目.

2.3 界面連續條件處理

圖3 為兩種不同力學參數飽和多孔介質的交界面在譜元離散情況下的相互作用力示意圖.在圖3中,介質交界面一條單元邊上有ngll個不等間距的計算節點(圖中圓點ngll=5);FRs和FRw分別為單元之間作用在固相和液相上的相互作用力;下標中的N和T分別表示交界面節點處作用力的法向分量和切向分量;A和B及對應的上標分別表示圖中上下介質交界面任意一個計算節點;e1,e2,e3,e4 為單元編號.

圖3 界面連續條件的譜元離散Fig.3 Spectral-element discretization of the continuous condition on the interface

結合式(6)交界面節點的應力和位移連續條件,可得到對應的相互作用力與位移連續條件,即

將上述連續條件與節點A和B的固液兩相動力平衡方程(參考式10)結合,即可求得圖3 中的8 個相互作用力.

需要說明的是,當圖3 中的飽和多孔介質向理想流體或彈性固體退化時,式(18)的界面連續條件亦會發生變化[1-2,5,36].將不同介質交界面的連續條件與對應節點的動力平衡方程相結合,即可求出對應的界面相互作用力.

2.4 人工邊界條件與地震動輸入

本文使用廖振鵬等[42-43]提出的多次透射公式(multi-transmitting formula,MTF)作為局部場地模型的人工邊界,該邊界具有簡單、通用、精度較高的特點.在飽和多孔介質的場地模型中,多次透射邊界需要同時應用于邊界節點固相位移和液相位移的x,y兩個方向分量之上,即

式中,N為透射階次,為二項式系數,以為例,其表示參考點j的x向固相位移在t=pΔt時刻的波場值,由下式定義

其中,Δt為時間步距,ca為人工波速.這里人工波速指的是MTF 邊界所使用的計算波速參數,理論上它可以設定為任意值,在實際計算中通常取為等于或略大于介質物理波速的某個值,本文取為固相縱波波速和液相縱波波速中的大值.

式(20)中參考點的位移通常不在波場的計算節點上,需要通過空間內插來獲得,而插值形式又與單元場函數的形式有關.采用本課題組開發的MTF 譜元離散格式實現[38-40],其特點為精度高且插值基函數與單元形函數一致.以為例,對應的MTF 譜元離散格式如下

表1 透射邊界(MTF)譜元離散格式中插值點的局部坐標Table 1 Local coordinates of the interpolation points in the spectral-element formulation of multi-transmitting formula

針對外源問題,地震動的輸入需要通過人工邊界來實現,這也是MTF 的一個優勢,應力型邊界或PML 邊界都不便于實現在人工邊界上的地震動輸入.MTF 作為位移型人工邊界,能方便地對全波場進行入射波場與散射波場的分離.圖4 為譜元模型中多次透射公式人工邊界條件和地震動輸入示意圖.其中,等號右邊的節點代表之前某時刻求得的全波場解;加號右邊節點代表已知的入射波場,本文通過傳遞矩陣方法(transfer matrix method,TMM)[41-44]求得;全波場解減去入射波場解,則得到之前某時刻的散射波場,即圖中加號左邊圓形實心節點和矩形實心節點,空心圓點表示MTF 的空間插值參考點,其值由式(21c)插值而得.求得MTF 各參考點的位移后,根據式(21a)可求得p+1 時刻的邊界節點位移,再結合波動方程算出的其余內節點位移,則得到p+1時刻完整的全波場解.

圖4 多次透射公式人工邊界條件與地震動輸入Fig.4 Artificial boundary condition and seismic wave input based upon multi-transmitting formula

3 算例驗證

選取了二維水平成層、不規則層狀界面和任意形狀界面的三種場地模型,介質皆為理想流體-飽和兩相土體-彈性固體,具體材料參數見表2.二維模型的側邊界和底面采用多次透射人工邊界(MTF,本文取為一階)來模擬無限域,上表面為自由地表,側邊界的入射波場采用傳遞矩陣方法(TMM)[41,44]得到的自由場,底面的入射波場則為圖5 所示的脈沖波.

圖5 輸入樣條脈沖波Fig.5 The input wave of a spline impulse

表2 介質參數表Table 2 Parameters of medium

3.1 水平層狀場地

圖6 為二維水平成層場地模型,下述所有工況的譜單元尺寸皆滿足 ΔxE＜λmin(λmin為三種介質中的最短波長),時間步距為 Δt=0.000 2 s,滿足穩定性條件.

圖6 二維水平成層場地模型Fig.6 Two-dimensiona model of horizontal layered site

本文方法(S E M) 除了與傳遞矩陣解析解(TMM)作對比外,還與有限元法(FEM)和統一框架計算理論相結合的結果作對比.其中,有限單元尺寸為 Δx＜λmin/10,對應的時間步距為Δt=0.000 05 s(Δx=1 m)或 Δt=0.000 2 s(Δx≥2 m).測點A,B,C,D位于底部邊界和介質交界面的中點處.

圖7 為P 波垂直入射情況下各測點的豎向位移.從圖7 中不難看出,本文方法與傳遞矩陣解析解在所有測點的位移都完全重合.有限元法與傳遞矩陣解析解在測點A和測點B的位移幾乎重合(見圖7(a)～圖7(d)),但在C點和D點的位移則存在一定誤差(見圖7(e)～圖7(h)),尤其是C點的飽和土固相位移和液相位移(見圖7(e)和圖7(f)).值得注意的是,譜元模型的計算節點總數為2499 個,而有限元模型的計算節點總數為14 883 個,這意味著譜元模型用比有限元模型少約83%的計算節點數,達到了比有限元模型更高的精度.

圖7 各測點的豎向位移Fig.7 Vertical displacement at measuring points

通過上述時程曲線觀察得出的譜元法模擬結果要優于有限元法模擬結果,其本質原因在于高階單元比低階單元能夠在更寬的頻率范圍內精確地模擬波場.數值離散網格對于波傳播過程而言相當于低通濾波器,高階方法的優勢在于它能夠達到的上限頻率更高.

為定量地考察這種頻率范圍的不同,這里以解析解作為參照,探討譜元法和有限元法數值解相對于解析解的譜比的差異.譜比計算公式如下

其中,S(f)為譜比,U(f)為頻譜幅值.

圖8 繪出了在不同網格尺寸條件下譜元法和有限元法計算出的C點固相、液相位移相對于解析解的譜比曲線.若譜比值越接近于1,表示該頻率處數值解越精確;反之,若偏離程度越大,則表示該頻率處數值解的誤差越大.圖中結果清楚地表明,譜元法能夠在比有限元法寬得多的頻帶范圍內給出精確的模擬結果.本算例中譜元法(四階單元)能夠精確模擬的波動頻率上限約為13 Hz,而有限元法(一階單元)只能達到約5 Hz;譜元法在0～13 Hz 范圍內模擬精度極好,而有限元法在0～5 Hz 范圍內仍存在少量誤差.此外,圖中結果還顯示出網格尺寸對波動模擬精度的影響規律,比較FEM 網格尺寸分別為1 m,2 m,4 m 的譜比曲線可知,網格尺寸越小,譜比誤差越低,即模擬精度越高.但是,縮小網格尺寸對模擬精度的提高遠不如提高單元階次所帶來的收益.總之,譜元法是適用于波動數值模擬的一種特殊的高階有限元,本算例表明高階單元的最大優勢在于能夠在寬得多的頻帶范圍內實現對波場的精確模擬.

圖8 不同方案情況下C 點的位移譜比Fig.8 Spectral ratio of displacement at point C with different scheme

需要說明的是,這里主要體現的是空間離散方式對波場頻散的影響,時步積分方法亦能產生較大影響[45-46].統籌考慮時空離散方法對數值模擬的頻散影響具有重要研究意義.鑒于本文的主題是使用顯式譜元模擬方法求解大型場地的地震動反應,這里不再深入討論.

3.2 不規則層狀界面場地

圖9 為不規則層狀界面場地模型,P 波從模型右下角以 θ=10°的傾斜角斜入射至波場中,介質分界面處的斜坡坡角為 α=14°.譜單元尺寸為ΔxE=2～20 m,時間步距為Δt=0.000 1 s.測點A,B,C,D位于介質交界面的角點處.

圖9 不規則層狀界面模型Fig.9 Irregular layered interface model

圖10 為各測點的位移時程圖.其中,時程曲線名稱對應的符號u和U分別表示固相位移和液相位移,上標b,s,w分別表示基巖、飽和土、理想流體,下標x,y分別表示測點的橫向與豎向.從圖10中可以看出,由于P 波斜入射以及存在不規則界面的原因,各測點的橫向位移出現了不同的波動.基巖和飽和土中的測點(A-D)由于處于坡度較小的坡面上,其橫向位移主要來源于入射波的水平分量,幅值較小,見圖10(a)～圖10(f).然而,在理想流體中,測點C與測點D的液相位移出現了較大的水平分量(見圖10(g)和圖10(h)).這主要是因為這兩個測點位于不規則界面的角點處,且理想流體的厚度相對較小,容易導致散射波的疊加,并在界面之間來回反射.

在圖10 中,測點C和測點D的位移曲線沒有觀察到一致性,主要是因為界面的位移連續性是針對界面的法向和切向的,而不是橫向與豎向.在實際編程過程中,界面測點的法向方向為界面上該點所在角的角平分線,切向方向則垂直于法向方向.下面將測點C和測點D的位移沿著法向方向進行分解.根據理想流體和飽和土的位移法向連續條件,有

圖10 各測點的位移Fig.10 Displacement at measuring points

其中,下標N表示測點的位移矢量沿著法向的分量,1 表示理想流體,2 表示飽和土,nx和ny表示單位法向矢量在橫向和豎向的分量.

同理,對于測點A和測點B,根據彈性固體和飽和土的位移法向連續條件,有

其中,2 表示飽和土,3 表示彈性固體(基巖).

圖11 為結合式(23)與式(24)計算得到的各測點法向位移.從圖11 中可以看出,彈性固體與飽和土分界面上的測點A和B以及飽和土與理想流體分界面上的測點C和D皆滿足位移法向連續條件.

圖11 各測點的法向位移Fig.11 Normal displacement at measuring points

圖12 給出了土體(基巖與飽和土)表面各計算節點固相位移峰值與輸入波幅值的比值,即|max(ut)/AP|.這里 max(ut) 為各節點的位移時程峰值,AP為輸入P波的幅值,圖中分別表示固相橫向位移與固相豎向位移.

圖12 土體表面的位移幅值Fig.12 Surface displacement amplitude of soil

從圖12 中可以看出,由于土體表面整體上比較平穩,且入射角較小,各節點的固相水平位移幅值較小,而豎向位移峰值較大;基巖表面的豎向位移幅值比飽和土表面的豎向位移幅值要大,這主要是由傳播介質的不同以及地形變化所導致的.

3.3 任意形狀界面場地

圖13 為任意形狀界面場地模型,P 波從模型右下角以 θ=5°的傾斜角斜入射至波場中.譜單元尺寸為 ΔxE=3～10 m,時間步距為 Δt=0.000 2 s.測點A,B,C,D位于介質交界面的拐點處.

圖13 任意形狀界面模型Fig.13 Arbitrary shape interface model

結合式(23)和式(24),圖14 給出了各測點的法向位移.從圖中可以看出,不同分界面上各測點的法向位移皆滿足位移連續條件.

圖14 各測點的法向位移Fig.14 Normal displacement at measuring points

圖15 給出了飽和土表面各計算節點固相位移峰值與輸入波幅值的比值.從圖中可以看出,由于飽和土表面起伏較大,部分節點的固相橫向位移的成分較大,接近0.5;固相豎向位移峰值隨節點的起伏分布有一定變化,但都在1.3 左右.

圖15 飽和多孔介質表面的位移幅值Fig.15 Surface displacement amplitude of saturated porous medium

3.2 節與3.3 節兩個模型的計算結果證明了本文方法應用于不規則界面復雜介質模型(理想流體-飽和多孔介質-彈性固體)的適用性,可用于大型河流湖泊、海洋區域工程、河谷盆地及大壩工程等場地的地震動數值模擬研究.

4 結論

如何對涉及到理想流體、飽和多孔介質、彈性固體三種介質的流固耦合地震波動問題進行高效的數值模擬是大型河流、湖泊以及海洋范圍內各種重要工程的抗震設防迫切需要解決的關鍵問題.本文基于飽和兩相介質波動的Biot 方程及其向理想流體和彈性固體的兩種退化形式方程,采用兼具波場模擬的高精度特性(源于譜插值)與復雜區域建模靈活性(源于與有限元相同的計算模式)的譜元法,結合簡單高效的廖氏透射邊界條件,建立了一種流固耦合地震波動問題的高效數值模擬方法.通過含有水平界面、不規則層狀界面和任意界面場地的數值模型算例,證明了本文方法的可行性,并且揭示出與目前大量使用的有限元法相比,該方法能夠在寬得多的頻帶范圍內實現對流固耦合場地中地震動的有效模擬.

本文方法除了具備現有流固耦合地震波動有限元模擬方法中完全顯式數值計算格式的一般特征之外,還具有從理論上嚴格地導出集中質量矩陣和高精度譜單元能夠很好地平衡地震動模擬精度與大規模數值計算效率這兩方面的優勢.文中討論主要針對二維問題進行,可以方便地推廣至三維情形.

此項研究對于數值模擬得到的地震動的寬頻帶特性分析具有理論意義,并且對于大型河流、湖泊以及海洋區域的場地地震反應分析具有工程應用價值.后續工作可以繼續研究大尺度譜單元對局部不均勻介質的有效建模、高精度人工邊界、地震動輸入方式以及時步積分方法等對地震動不同頻段特性的影響,或將各向異性介質[15,47-48]納入統一框架理論,更重要的是考慮飽和兩相介質的非線性并使用并行計算技術開發實用的三維計算程序,為解決實際工程問題服務.