李偉華， 王文強(qiáng)

(1.北京交通大學(xué) 土木建筑工程學(xué)院，100044；2. GeoApplication Engineers Pte Ltd, 367903, Singapore)

動荷載下地基振動響應(yīng)與波動規(guī)律作為土動力學(xué)基本問題之一，很早就被學(xué)者關(guān)注[1]。近年來隨著我國高速鐵路的發(fā)展，為這一問題的研究注入了新的活力，并促使其發(fā)展和完善。最初關(guān)于地基振動響應(yīng)問題的研究是在單相彈性半空間中展開，并取得了相當(dāng)豐碩的成果。自從Biot[2]建立了飽和多孔介質(zhì)的波動方程后飽和土地基的動力問題研究得到了長足的發(fā)展。Philippacopoulos[3]、Rajapakse[4]、陳龍珠[5]、金波[6]、陳勝立[7]等就此發(fā)表了一系列論文，逐步完善了地基與基礎(chǔ)的動力相互作用問題。然而實(shí)際上地球表面的土體大多數(shù)處于非飽和狀態(tài)，通常我們所說的飽和土，實(shí)際上會或多或少含有一些微量氣泡，是一種準(zhǔn)飽和土。已有的研究證明，即使微量氣泡的存在也會對土中波的傳播特性產(chǎn)生很大的影響[8]，所以研究非飽和土地基的動力響應(yīng)更具有實(shí)際意義。

非飽和土體的動力特征比飽和土體動力特征要復(fù)雜的多，現(xiàn)有研究多停留在非飽和多孔介質(zhì)中波的傳播特性上。Berryman[9]，Santos[10]，Lo[11]，Tuncay[12]，Vardoulakis[13]，Wei[14]和Lu[15]等分別依據(jù)Biot理論、混合物理論或二者的結(jié)合各自建立了非飽和多孔介質(zhì)的動力方程，確定非飽和多孔介質(zhì)中存在著1種剪切波和3種壓縮波(P1，P2和P3波)，并研究了其傳播特征。由于非飽和介質(zhì)動力方程的復(fù)雜性，往往需要借助數(shù)值方法對方程進(jìn)行了求解。Zhao等[16]采用交叉網(wǎng)格有限差分方法求解了Santos提出的非飽和多孔介質(zhì)方程，并結(jié)合完美匹配層邊界研究了非飽和多孔介質(zhì)中的波動。Wei[17]和Ravichandran[18]根據(jù)Wei所建立方程，建立有限元方法并分析非飽和土壩的地震響應(yīng)。Gatmiri[19]和Maghoul[20]在準(zhǔn)靜態(tài)條件下推導(dǎo)了非飽和多孔介質(zhì)動力方程的邊界積分表達(dá)式，并給出方程的2D和3D的時域、頻域解。Moldovan等[21]采用Hybrid-Trefftz應(yīng)力單元法求解了Wei所建立的方程。目前，關(guān)于非飽和土地基的動力響應(yīng)問題的研究并不多見。徐明江[22]以三相多孔介質(zhì)模型為基礎(chǔ)建立非飽和土的動力控制方程，采用解析方法研究了非飽和土半空間在表面荷載作用下的穩(wěn)態(tài)和瞬態(tài)動力響應(yīng)，以及非飽和土地基上剛性和彈性圓板基礎(chǔ)的垂直振動問題。

鑒于用解析法求解非飽和土地基的動力響應(yīng)問題的諸多限制，只有在少數(shù)特定的邊界條件下才能求得，本文建立了一種非飽和多孔介質(zhì)動力分析顯式有限元方法。首先在Wei所建立的非飽和多孔介質(zhì)動力方程的基礎(chǔ)上，將其中的土水特征曲線模型改為能夠更好描述接近飽和時土體特征的van Genuchten模型，然后通過時空離散，建立該方程的顯式有限元表達(dá)式。該方法采用解耦計(jì)算，無需求解聯(lián)立方程組，能夠有效節(jié)省內(nèi)存空間和提高計(jì)算速度，在求解Biot飽和多孔介質(zhì)波動方程中已經(jīng)得到了很好的應(yīng)用[23]。利用此方法與局部透射人工邊界結(jié)合，編制計(jì)算程序，分別建立剛性和柔性條形基礎(chǔ)作用下非飽和土地基模型，分析了在豎向均布脈沖荷載作用下，飽和度、孔隙率等因素對的非飽和土地基振動響應(yīng)的影響。

1 非飽和土地基動力方程

Wei以多孔介質(zhì)混合理論為基礎(chǔ)，考慮動態(tài)相容條件建立非飽和多孔介質(zhì)波動方程：

(1a)

(1b)

(1c)

式中：上角s，w和a分別代表非飽和土的固相、液相和氣相。nα，ρα，uα(α=s,w,a)分別表示α相的體積分?jǐn)?shù)、密度以及位移。rαα為彈性系數(shù)，文獻(xiàn)[24]給出其與非飽和土各相的體積模量Kα、固體骨架的體積模量K、剪切模量G及有效應(yīng)力系數(shù)αB的關(guān)系，

(2a)

(2b)

(2c)

(2d)

(2e)

(2f)

其中，

(3a)

(3b)

Θw和Θa是描述毛細(xì)管平衡條件的參數(shù)，對給定的飽和度Sr，有：

(4)

式中：pa-pw為基質(zhì)吸力，對于給定飽和度，由土水特征曲線模型確定。Wei建立的非飽和多孔介質(zhì)波動方程中，土水特征曲線及非飽和土的滲透系數(shù)均采用Brooks-Corey模型[25]。該模型存在導(dǎo)數(shù)不連續(xù)點(diǎn)，不能很好地描述土體接近飽和時的特征，而接近飽和時土水特征曲線形狀上的微小變化都會導(dǎo)致預(yù)測水力傳導(dǎo)率函數(shù)的顯著變化，從而影響非飽和帶水流模擬的結(jié)果，同時也會造成數(shù)值計(jì)算困難。為了更好地描述土體接近飽和時的特征，本文將其中的土水特征曲線及非飽和土的滲透系數(shù)改為van Genuchten提出的經(jīng)驗(yàn)公式[26]。則基質(zhì)吸力表示為，

(5)

式中：α，m，d為模型參數(shù)，且d=1-1/m，

(6)

式中：SrW為束縛飽和度；SrN為空氣殘余飽和度。

由線彈性狀態(tài)下非飽和多孔介質(zhì)的應(yīng)力-應(yīng)變關(guān)系，可以得到Θw和ΘN的關(guān)系，

(7)

聯(lián)立式(4)～(7)即可求得Θw和Θa。

參數(shù)ξf(f=w,a)與滲透率k的關(guān)系如下，

(8)

(9)

2 顯式有限元方程的建立

2.1 非飽和多孔介質(zhì)動力方程的矩陣形式

定義如下矩陣：

(10)

(11)

(12)

利用上述矩陣表達(dá)，方程(1)可以表示為：

(13)

其中

2.2 有限元方程列式的建立

設(shè)N為單元形函數(shù)矩陣，則單元的位移向量uα,e和力向量fα,e可以分別表示為：

(14)

(15a)

(15b)

(15c)

式中：應(yīng)變矩陣B=L·N。

2.3 離散方程及解耦方法

對于二維問題，采用4節(jié)點(diǎn)4邊形等參單元，節(jié)點(diǎn)i和與其直接相鄰的節(jié)點(diǎn)k構(gòu)成局部節(jié)點(diǎn)系，此節(jié)點(diǎn)系中包含9個節(jié)點(diǎn)。根據(jù)式(13)可以直接得到如下離散公式：

(16a)

(16b)

(16c)

其中，

(17)

(18)

(19)

(20)

(21)

2. 4 節(jié)點(diǎn)動力反應(yīng)的顯式表達(dá)式

利用有限元方法建立離散的動力方程后，采用中心差分法與Newmark平均加速度近似格式[23]相結(jié)合的方法求解節(jié)點(diǎn)動力方程式(16)。得到由前一時刻的位移和速度表示的當(dāng)前時刻位移的顯式表達(dá)式，

(22a)

(22b)

(22c)

由前一時刻的位移和速度表示的當(dāng)前時刻的節(jié)點(diǎn)速度反應(yīng)的顯式表達(dá)：

(23a)

(23b)

(23c)

由前一時刻的速度和此時刻的速度表示的加速度的顯式表達(dá)式

(24)

式(22)～(24)為非飽和多孔介質(zhì)局部節(jié)點(diǎn)系下某一節(jié)點(diǎn)的動力反應(yīng)的顯式表達(dá)式，該節(jié)點(diǎn)的動力反應(yīng)只與相鄰節(jié)點(diǎn)的動力反應(yīng)有關(guān)，從而實(shí)現(xiàn)了方程的解耦求解。

2.5 動力人工邊界

對半無限域地基進(jìn)行有限元模擬，需要引入動力人工邊界，本文采用局部透射邊界[27]，非飽和多孔介質(zhì)中的4種波(P1，P2，P3，S波)都以其固有的波速向前傳播。無論是固相、液相還是氣相，對某一特定的波來說，都是以同一波速沿同一方向向前傳播，且氣相、液相位移與固相位移具有固定的關(guān)系。因此，可以假設(shè)非飽和多孔介質(zhì)的固相、液相、氣相位移分別以同一波速通過人工邊界，針對固相、液相、氣相分別應(yīng)用多次透射人工邊界。

3 有限元求解精度驗(yàn)證

鑒于目前關(guān)于非飽和多孔介質(zhì)動力反應(yīng)問題分析的解析解較少，且非飽和多孔介質(zhì)動力反應(yīng)的控制方程不統(tǒng)一，很難找到與本文中的控制方程符合的解析解，因而為了驗(yàn)證計(jì)算方法的可行性，本文采用了以下兩種方法：方法一：將非飽和多孔介質(zhì)動力控制方程退化為飽和多孔介質(zhì)，與文獻(xiàn)[28]的解析解進(jìn)行對比驗(yàn)證。假設(shè)孔隙中只含有水，即非飽和多孔介質(zhì)的飽和度為1(氣相體積分?jǐn)?shù)為0)，將氣相相關(guān)參數(shù)設(shè)為零，方程(1)可以退化到飽和狀態(tài)：

(25a)

(25b)

式中的各項(xiàng)參數(shù)與Biot提出的飽和多孔介質(zhì)動力方程的參數(shù)對比情況見表1。

表1 參數(shù)對比

表2 飽和土材料參數(shù)

圖1 驗(yàn)證模型及網(wǎng)格劃分示意圖

圖2 驗(yàn)證一： 固相位移時程圖

圖3 驗(yàn)證一：液相位移時程

方法二：將氣相用液相代換，且假設(shè)兩種液體的材料參數(shù)相同，這樣可以保證方程的完整性。采用和方法一相同的分析模型和參數(shù)，圖4給出了孔隙被兩種相同液體充滿時，自由表面固相位移隨時間的變化曲線。可以看出，第二種方法的計(jì)算結(jié)果與文獻(xiàn)[28]的解析解十分接近，由于將單種液體劃分為兩種液體，計(jì)算過程中多考慮一個界面的影響，會使計(jì)算結(jié)果產(chǎn)生一定差別。但結(jié)果依舊能很好地反映出本文計(jì)算方法的可靠性。

圖4 驗(yàn)證二：固相位移時程圖

4 非飽和土地基振動分析

4.1 地基模型及土層參數(shù)

在地基表面處有一寬為4 m(2b=4 m)的條形基礎(chǔ)，基礎(chǔ)底面作用幅值σ0=1 kPa作用時間T=0.5 s的豎向均布脈沖荷載。分別考慮剛性基礎(chǔ)和柔性基礎(chǔ)兩種情況，不考慮地基與基礎(chǔ)的相互作用及基礎(chǔ)自重，將基礎(chǔ)簡化為豎向均布條形荷載作用在地基表面。計(jì)算模型如圖5所示。取計(jì)算區(qū)域?yàn)?0 m×30 m，網(wǎng)格尺寸為1 m×1 m，網(wǎng)格劃分及坐標(biāo)如圖所示。地基土的物理性質(zhì)見表3。在基礎(chǔ)中心點(diǎn)A以下不同深度及地基表面距A點(diǎn)不同距離處取觀察點(diǎn)，分析非飽和土地基振動規(guī)律及飽和度對非飽和土地基振動的影響。

圖5 地基計(jì)算模型圖

K/PaG/Paρs/(kg·m-3)ρw/(kg·m-3)ρa(bǔ)/(kg·m-3)2.59×1071.94×107270010001.10Ks/PaKw/PaKa/Pavw/(Pa·s)va/(Pa·s)3.6×10102.0×1091.0×1051.0×10-31.5×10-5k/μm2αBαmn01.01.01.0×10-42.00.36或0.6

4.2 剛性基礎(chǔ)作用下非飽和土地基動力反應(yīng)分析

考慮非飽和土地基表面與基礎(chǔ)接觸面上各點(diǎn)滿足以下條件：①剛性接觸面上各點(diǎn)固相豎向位移相等；②剛性接觸面內(nèi)外均為透水邊界。

圖6給出了當(dāng)土體的初始孔隙率n0=0.36時，基礎(chǔ)中心點(diǎn)下各點(diǎn)的固相豎向位移幅值，其中圖6(a)為各點(diǎn)豎向位移隨飽和度變化的關(guān)系曲線；圖6(b)為不同飽和度下基礎(chǔ)中心點(diǎn)下豎向位移沿深度(z/b)變化曲線。

圖6 剛性基礎(chǔ)中心點(diǎn)下不同深度處各點(diǎn)豎向位移幅值(n0=0.36)

從圖6中可以看出，當(dāng)?shù)鼗敛伙柡?Sr<1)時，隨著飽和度的增加，剛性基礎(chǔ)中心點(diǎn)下各點(diǎn)豎向位移減小，且減小程度隨著飽和度的增加變得明顯。這主要是由于隨著飽和度增加，土中基質(zhì)吸力降低，有效壓應(yīng)力減小所致。而且當(dāng)土體飽和度較小時，孔隙中含有較多的空氣，孔隙流體的壓縮性增加，由荷載引起的孔隙流體壓力較小，飽和度變化對骨架位移影響不大。當(dāng)?shù)鼗镣耆柡?Sr=1)時，各點(diǎn)豎向位移較Sr<1時有顯著降低。說明在土體接近完全飽和時，即使飽和度發(fā)生微小變化都會影響基礎(chǔ)中心點(diǎn)下各點(diǎn)的豎向振動。這是因?yàn)楫?dāng)土體接近飽和時，土中氣體以封閉氣泡形式存在，氣體含量的減小，會導(dǎo)致孔隙流體的體積模量顯著增加，從而使得土體抵抗變形能力增加。從圖中還可以看到，基礎(chǔ)中心點(diǎn)下各點(diǎn)的位移沿深度方向衰減?？拷行狞c(diǎn)處衰減較快，遠(yuǎn)離中心點(diǎn)各點(diǎn)位移隨深度變化趨于平緩。

圖7給出了當(dāng)土體的初始孔隙率n0=0.36時，地表固相豎向位移幅值，其中圖7(a)為地表各觀察點(diǎn)豎向位移隨飽和度變化的關(guān)系曲線；圖7(b)為不同飽和度下地表豎向位移沿水平方向(x/b)變化曲線。

圖7 剛性基礎(chǔ)地基表面各觀察點(diǎn)豎向位移幅值(n0=0.36)

從圖7中可以看出，和基礎(chǔ)中心點(diǎn)下豎向位移一樣，隨著地基土體飽和度的增加，地表各點(diǎn)豎向位移減小。當(dāng)飽和度較小時，各點(diǎn)豎向位移隨飽和度變化不大，飽和度較大時，飽和度對各點(diǎn)豎向位移的影響變得愈加明顯。且當(dāng)Sr=1時，各點(diǎn)豎向位移較Sr<1時有顯著降低。從圖中還可以看到，地表各觀察點(diǎn)的豎向位移隨觀察點(diǎn)到荷載中心點(diǎn)距離的增加而減小。因?yàn)槭莿傂曰A(chǔ)，荷載作用范圍內(nèi)各點(diǎn)位移相等，當(dāng)12時變化逐漸平緩。

圖8給出了當(dāng)土體的初始孔隙率n0=0.36時，地表固相水平位移幅值，其中圖8(a)為地表各觀察點(diǎn)水平位移隨飽和度變化的關(guān)系曲線；圖8(b)為不同飽和度下地表水平位移沿水平方向(x/b)變化曲線。圖中負(fù)號表示水平位移與坐標(biāo)方向相反。

圖8 剛性基礎(chǔ)地基表面各觀察點(diǎn)水平位移幅值(n0=0.36)

從圖8中可以看出，當(dāng)Sr<1時，隨著飽和度的增加，地表各點(diǎn)水平位移大小增加，飽和度較小時，變化不大，飽和度較大時，增加明顯。這主要是因?yàn)楫?dāng)Sr<1時，隨著飽和度的增加，土體中的基質(zhì)吸力降低，水平方向有效拉應(yīng)力增加，導(dǎo)致水平方向產(chǎn)生的與坐標(biāo)方向相反的位移大小增加。當(dāng)Sr=1時，地表各點(diǎn)水平位移大小較Sr<1時反而降低。這是由于孔隙中封閉的氣泡使得地基土體接近飽和時孔隙流體的體積模量較飽和時顯著降低。從圖中還可以看到，地表同一觀察點(diǎn)的水平位移較豎向位移小近一個數(shù)量級，地表各觀察點(diǎn)的水平位移隨觀察點(diǎn)到荷載中心點(diǎn)距離的增加先增加后減小，在x/b=2時達(dá)到最大值。

為了分析土體初始孔隙率不同時飽和度對地基振動的影響，圖9～圖11分別給出了當(dāng)土體初始孔隙率n0=0.6，其他參數(shù)保持不變時，不同飽和度下，基礎(chǔ)中心點(diǎn)下以及地表的固相豎向、水平位移幅值分布。從圖中可以看出，n0=0.6時飽和度對各點(diǎn)振動響應(yīng)的影響規(guī)律與n0=0.36時相同。當(dāng)土體其他參數(shù)相同時，初始孔隙率大的，飽和度對各點(diǎn)振動響應(yīng)的影響較大。當(dāng)Sr<1時，初始孔隙率大，地基土體各點(diǎn)振動反應(yīng)較初始孔隙率小時增加明顯；當(dāng)Sr=1時，初始孔隙率的大小對地基振動響應(yīng)影響不大。再次說明非飽和土中氣相狀態(tài)的存在對地基振動響應(yīng)的影響。

圖9 剛性基礎(chǔ)中心點(diǎn)下豎向位移幅值沿深度分布曲線(n0=0.6)

圖10 剛性基礎(chǔ)地基表面豎向位移沿水平方向分布曲線(n0=0.6)

圖11 剛性基礎(chǔ)地基表面水平位移沿水平方向分布曲線(n0=0.6)

4.3 柔性基礎(chǔ)作用下非飽和土地基動力反應(yīng)分析

考慮非飽和多孔介質(zhì)的自由面與柔性基礎(chǔ)接觸面上各點(diǎn)滿足以下條件：① 柔性接觸面上各點(diǎn)的位移不加以約束；②考慮柔性接觸面間的滑動，在接觸面內(nèi)外皆為透水邊界。

圖12給出了不同飽和度下，土體初始孔隙率分別為0.36、0.6時，基礎(chǔ)中心點(diǎn)下固相豎向位移幅值沿深度的分布。

圖12 柔性基礎(chǔ)中心點(diǎn)下豎向位移幅值沿深度分布曲線

從圖12中可以看出，與剛性基礎(chǔ)相比，柔性基礎(chǔ)條件下，基礎(chǔ)中心點(diǎn)下相同位置處的豎向位移幅值要小。這與現(xiàn)有的研究結(jié)論，隨著基礎(chǔ)柔度增加，地基動力柔度系數(shù)增加是一致的[29]。在柔性基礎(chǔ)條件下，當(dāng)Sr<1時，飽和度的變化對基礎(chǔ)中心點(diǎn)下豎向位移幅值影響不如剛性基礎(chǔ)條件下明顯，初始孔隙率較大時，飽和度對基礎(chǔ)中心點(diǎn)下固相位移幅值影響稍有增加，隨著飽和度增加基礎(chǔ)中心點(diǎn)下固相位移幅值減小。與剛性基礎(chǔ)相同的是，Sr=1時各點(diǎn)豎向位移較Sr<1時有顯著降低；當(dāng)Sr<1時，初始孔隙率大，地基土體各點(diǎn)振動反應(yīng)較初始孔隙率小時增加明顯；當(dāng)Sr=1時，初始孔隙率的大小對地基振動響應(yīng)影響不大。

圖13給出了不同飽和度下，土體初始孔隙率分別為0.36、0.6時，地基表面固相豎向位移幅值沿水平方向的分布曲線。

從圖13中可以看出，與基礎(chǔ)中心點(diǎn)下豎向位移類似，柔性基礎(chǔ)地基表面的豎向位移幅值較剛性基礎(chǔ)條件下相同位置處的要小；當(dāng)Sr<1時，飽和度的變化對地表豎向位移幅值影響不如剛性基礎(chǔ)條件下明顯，初始孔隙率較大時，飽和度影響稍有增加，隨著飽和度增加地表同一點(diǎn)處的豎向位移幅值減小；與剛性基礎(chǔ)相同的是，Sr=1時各點(diǎn)豎向位移較Sr<1時有顯著降低；當(dāng)Sr<1時，初始孔隙率大，地基土體各點(diǎn)振動反應(yīng)較初始孔隙率小時有所增加；當(dāng)Sr=1時，初始孔隙率的大小對地基振動響應(yīng)影響不大。因?yàn)槭侨嵝曰A(chǔ)，荷載作用范圍內(nèi)地表各點(diǎn)豎向位移隨x/b的增加而衰減，并在0.5

圖13 柔性基礎(chǔ)地基表面豎向位移幅值沿水平方向分布曲線

圖14 柔性基礎(chǔ)地基表面水平位移幅值沿水平方向分布曲線

圖14給出了不同飽和度下，土體初始孔隙率分別為0.36、0.6時，地基表面固相水平位移幅值沿水平方向的分布曲線。從圖中可以看出，柔性基礎(chǔ)地基表面的水平位移幅值較剛性基礎(chǔ)條件下相同位置處的??；當(dāng)Sr<1時，在兩種初始孔隙率條件下，飽和度的變化對地表水平位移幅值影響均不明顯。Sr=1時各點(diǎn)水平位移較Sr<1時有顯著降低，尤其是在荷載作用范圍外(x/b>1)；初始孔隙率不同，地基表面相同點(diǎn)處水平位移幅值差別不大。說明柔性基礎(chǔ)條件下，飽和度及初始孔隙率對地表水平位移影響有限。

5 結(jié) 論

針對目前非飽和土地基的振動響應(yīng)研究剛剛起步，分析方法尚不完善，本文在現(xiàn)有的非飽和孔隙介質(zhì)波動方程的基礎(chǔ)上，采用解耦技術(shù)，建立非飽和多孔介質(zhì)動力反應(yīng)分析的顯式有限元方法。

采用該方法，分別建立剛性和柔性條形基礎(chǔ)作用下非飽和土地基振動模型，分析了在豎向均布脈沖荷載作用下，飽和度、孔隙率等因素對的非飽和土地基動力反應(yīng)的影響。計(jì)算表明：①剛性基礎(chǔ)條件下地基的振動明顯高于相同條件下柔性基礎(chǔ)地基振動，這與現(xiàn)有研究結(jié)果，隨著基礎(chǔ)柔性增加，地基動力柔度系數(shù)增加是一致的；②剛性基礎(chǔ)條件下，Sr<1時，隨著飽和度增加各點(diǎn)豎向位移減小，水平位移增大；Sr=1時各點(diǎn)豎向、水平位移均比Sr<1時明顯減小，說明非飽和土中氣相狀態(tài)的存在與否對剛性基礎(chǔ)條件下地基振動具有較為顯著的影響。③柔性基礎(chǔ)條件下，Sr<1時，飽和度對各點(diǎn)的豎向、水平位移影響有限；Sr=1時各點(diǎn)豎向、水平位移均比Sr<1時明顯減小。④飽和度對地基振動的影響與土體的孔隙率有關(guān)，其他條件相同時，初始孔隙率高的，飽和度影響較大，且相同條件下地基動力反應(yīng)大。

由于篇幅限制，本文只討論了飽和度對各點(diǎn)振動位移的影響，對于非飽和土地基其他振動響應(yīng)如動應(yīng)力及孔隙流體壓力分布等將在后續(xù)文章中討論。