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

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

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

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

引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

兩種不同力學參數(shù)飽和多孔介質(zhì)交界面的連續(xù)條件為:

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

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

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

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

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

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

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

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

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

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

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

2.1 譜元離散方程

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

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

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

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

形函數(shù)表達式為

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

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

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

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

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

2.2 時域逐步積分

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

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

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

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

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

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

其中,nc為與節(jié)點k相關(guān)的節(jié)點的數(shù)目.

2.3 界面連續(xù)條件處理

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

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

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

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

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

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

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

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

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

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

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

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

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

3 算例驗證

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

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

表2 介質(zhì)參數(shù)表Table 2 Parameters of medium

3.1 水平層狀場地

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

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

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

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

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

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

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

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

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

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

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

3.2 不規(guī)則層狀界面場地

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

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

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

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

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

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

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

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

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

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

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

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

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

3.3 任意形狀界面場地

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

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

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

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

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

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

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

4 結(jié)論

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

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

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