凌道盛 ，王 筠 ，單振東，丁皓江

（1. 浙江大學(xué) 軟弱土與環(huán)境土工教育部重點實驗室，杭州 310058；2. 浙江大學(xué) 巖土工程研究所，杭州 310058；3. 中國地震局工程力學(xué)研究所，哈爾濱 150080）

1 引 言

多孔介質(zhì)瞬態(tài)響應(yīng)問題的研究是巖土工程、地震工程和地球物理等學(xué)科領(lǐng)域的一個重要課題。在巖土工程中，如果多孔介質(zhì)內(nèi)充滿單一流體，則稱這種多孔介質(zhì)是飽和的；如果多孔介質(zhì)內(nèi)存在多種流體，則稱這種多孔介質(zhì)是非飽和的[1]。大多數(shù)情況下，多孔介質(zhì)內(nèi)有兩種流體共存。它們可能是兩種液體，如石油工程中的油和水；也可能是一種液體和一種氣體，如非飽和土中的水和氣。

自從Biot[2－3]1956年建立飽和多孔介質(zhì)波動理論以來，多孔介質(zhì)的波動理論引起了包括巖土工程在內(nèi)的多個領(lǐng)域的學(xué)者的廣泛重視。由于存在慣性、黏性和機械耦合作用，多孔介質(zhì)的動力問題比單相彈性介質(zhì)復(fù)雜得多，大多數(shù)問題只能通過時間和空間上離散等數(shù)值方法進行求解，即使是一維的瞬態(tài)響應(yīng)問題，也只能在某些特殊情形下才能給出解析解。Garg等[4]運用Laplace變換法求解一維瞬態(tài)響應(yīng)問題，得到了黏性耦合作用為0或無窮大兩種極端情況下時域內(nèi)的解析解；Simon等[5]運用Laplace變換給出了材料參數(shù)滿足動力相容條件時的解析解；Gajo等[6]給出了分離變量型的位移通解和飽和多孔介質(zhì)層上下表面同時作用階躍位移的精確解；Schanz等[7]運用 Laplace變換給出了黏性耦合為 0時的解析解；凌道盛等[8]給出了單層不可壓縮飽和多孔介質(zhì)一維瞬態(tài)響應(yīng)的精確解；Shan等[9]運用分離變量法給出了3類非齊次邊界條件的精確解；凌道盛等[10]提出了一種適用于單層飽和多孔介質(zhì)任意非齊次邊界條件的半解析方法。

針對非飽和多孔介質(zhì)一維瞬態(tài)響應(yīng)問題的研究成果相對較少，Li等[11]基于混合物理論，運用Laplace變換研究了單層問題，給出了變換域內(nèi)的解析解和時域內(nèi)的數(shù)值Laplace逆變換解；Shan等[12]給出了單層問題3類典型非混合邊界條件的級數(shù)解和半無限問題兩類邊界條件的積分解。1964年Brutsaert[13]指出，非飽和多孔介質(zhì)內(nèi)存在3種類型的壓縮波，這與飽和多孔介質(zhì)（兩類壓縮波）存在顯著的不同，深入認識非飽和多孔介質(zhì)一維瞬態(tài)響應(yīng)是理解多孔介質(zhì)內(nèi)波動傳播的基礎(chǔ)，有必要對非飽和多孔介質(zhì)一維瞬態(tài)響應(yīng)問題開展進一步的研究。

綜上可見，針對非飽和多孔介質(zhì)的一維瞬態(tài)響應(yīng)問題的研究還比較少，目前還沒有適用混合邊界的解析解，而且基于 Laplace變換法的求解方法存在數(shù)值 Laplace逆變換穩(wěn)定性較差、在響應(yīng)時間較長時可能出現(xiàn)發(fā)散等不足[14]。本文以一類混合邊界條件為例，給出了一種考慮孔隙流體和固體顆粒壓縮性以及慣性、黏滯和機械耦合作用的半解析解法。該解法不受邊界條件和材料參數(shù)的限制，具有很好的適用性。

2 基本方程和定解條件

本文研究的一維單層非飽和多孔介質(zhì)及其坐標系如圖1所示。圖中，H為非飽和多孔介質(zhì)層厚度。

2.1 基本方程

Zienkiewicz 等[15]和 Li等[16]將 Biot理論[2－3]推廣到非飽和多孔介質(zhì)領(lǐng)域，提出了非飽和多孔介質(zhì)波動理論的基本方程。對于多孔介質(zhì)孔隙內(nèi)部含有兩種不同流體，即濕相流體和干相流體的一維問題，方程可以表述為如下形式：

圖1 單層非飽和多孔介質(zhì)示意圖Fig.1 Sketch of single layer unsaturated porous media

式中：帶有下標w或n的變量分別是與濕相或干相流體有關(guān)的變量。σz為總應(yīng)力，以拉為正；Pw和Pn為濕相和干相流體的孔隙流體壓力，以壓為正，且有稱為有效孔隙流體壓力；Uz是固體骨架的絕對位移；Wz、Vz分別是濕相和干相的流體相對位移，且有和Vz=為流體的絕對位移；n為孔隙率；Sw和Sn分別為濕相流體和干相流體的飽和度，且有 Sw+Sn=1；ρs、ρw、ρn分別為固體顆粒、濕相流體和干相流體的密度，非飽和多孔介質(zhì)總體密度為分別為濕相流體和干相流體的動力滲透系數(shù)；γ和μ為固體骨架的 Lame常數(shù)；α、Qw和Qn均為與固體顆粒和流體壓縮性有關(guān)的系數(shù)，且有(α - n)/Ks+ n/Kn。其中Kb、Ks、Kw和Kn分別為固體骨架、固體顆粒、濕相流體和干相流體的體積模量，且有 Kb=λ+2μ / 3；t為時間。

式（1）、（2）、（3）分別為非飽和多孔介質(zhì)整體平衡方程和兩種流體的運動平衡方程；式（4）、（5）分別為兩種流體的滲流連續(xù)方程；式（6）為固體骨架的本構(gòu)方程；式（7）為固體骨架的幾何方程。

為方便起見，引入如下無量綱量：

利用上述無量綱量，基本方程（1）～（7）可以寫成如下矩陣形式：

其中：

其中：

將式（10）代入式（9）中，得到用位移矢量u表示的控制方程為

2.2 定解條件

限于篇幅，本文僅以如下混合邊界條件為例給出非飽和多孔介質(zhì)一維瞬態(tài)響應(yīng)半解析解。

式中： fij(τ)(i = 0,1; j = 1, 2,3)為可以具有任意形式的指定函數(shù)。

初始條件為

式中： gij(ζ)(i = 0,1; j = 1,2,3)為可以具有任意形式的指定函數(shù)。

3 問題的求解

3.1 邊界條件齊次化

將位移表達式分解為兩部分：

其中：

選取適當?shù)?us，使其滿足非齊次邊界條件（16）。本文選取的 us表達式為

其中：

將式（18）代入到式（15），得到 ud表示的控制方程、邊界條件和初始條件：

當ζ=0時，有

當ζ=1時，有

當τ=0時，有

3.2 特征值和特征函數(shù)

為了求解 ud，先求 ud對應(yīng)的無阻尼特征值問題，令

式中：ω≥ 0。式（24）～（26）對應(yīng)的無阻尼系統(tǒng)特征方程為

當ω=0時，式（29）具有如下形式的解：

將式（32）代入到式（30）、（31），可求得

式中：c1、c2為任意常數(shù)?？梢姦?0對應(yīng)的兩個模態(tài)分別為

當ω＞0時，式（29）具有如下形式的解：

式中：βd為待定常數(shù)，將式（35）代入式（29）得

為使方程（36）有非零解，其系數(shù)矩陣行列式必須為0。由此可得到如下特征方程：

其中：

由式（37）可求得6個特征值，記為βdi(i=1, 2,3, 4, 5, 6)，在一般情況下，βdi均為虛數(shù)，且有

式中：β1、β2、β3是只與材料參數(shù)有關(guān)，與ω值無關(guān)的常數(shù)。利用式（35），可求得特征值βdi對應(yīng)的特征向量，記為Φi(i=1, 2, 3, 4, 5, 6)，則有

其中：

此時X的表達式可以寫成：

將式（42）代入邊界條件（30）得

因為方程（43）的系數(shù)矩陣行列式不為 0，所以式（43）只有零解，即 a4=a5=a6= 0。此時，式（42）化簡為

將式（44）代入邊界條件（31）得

式中： i= 1, 2, 3。

若式（45）有非零解，其系數(shù)矩陣行列式必為零，即

上式即為確定ω的特征方程，由此可得到一系列的特征值，從小到大排列，記第 p(p ≥ 3)個特征值為ωp，記ωp對應(yīng)的特征向量為由式（45）可求得

因此，Xp的表達式就可以寫成：

此時所有ω對應(yīng)的模態(tài)都已經(jīng)求出。為方便起見，記 ω1=ω2= 0，令

則所有模態(tài)均可用式（49）表示，只是p=1,2,3…。

3.3 特征函數(shù)的正交性

設(shè)ωp和ωq為兩個不同的特征值，它們對應(yīng)的特征函數(shù)為Xp和Xq。利用M矩陣的對稱性質(zhì)和式（29）可得

采用分部積分技術(shù)并應(yīng)用邊界條件，容易證明上式第一個積分值為0。

將Xp和Xq代入可得

采用分部積分技術(shù)，有

根據(jù)邊界條件（30）、（31），上式可簡化為

用同樣的方法可以得到

于是，特征函數(shù)的正交性就得到證明。

當 ωp≠ωq時，有

特征函數(shù)M的正交性得到證明，為方便起見，將特征函數(shù)關(guān)于矩陣M規(guī)范化處理，定義規(guī)范化的特征矢量為

3.4 動力問題的解

對有阻尼系統(tǒng)的受迫振動，根據(jù)變異系數(shù)法，設(shè) ud具有以下形式的解：

代入式（24），并利用式（29）得

其中：

將式（64）代入初始條件（27），得

求解由式（66）、（70）構(gòu)成的初值問題，即可得到 ud，利用式（18）即可得到原問題u的解。

4 常微分方程組的精細時程積分解

常微分方程組（66）中含有無限個常系數(shù)微分方程，由于ω已經(jīng)按照從小到大的順序排列，因此，在實際求解時，可以根據(jù)荷載性質(zhì)及精度要求截取前N項進行計算。將式（66）改寫成如下矩陣形式：

對應(yīng)的初始條件仍為式（70）。

當濕相流體和干相流體的滲透系數(shù)都趨向無窮大時，cpq趨近于0，阻尼陣為零矩陣，式（71）解耦，此時可求得解析解。一般情況下，為滿陣，式（70）、（71）構(gòu)成的初值問題難以獲得精確解，本文采用鐘萬勰[17]提出的精細時程積分法求解上述初值問題。

5 數(shù)值算例

以下算例均假定多孔介質(zhì)內(nèi)部孔隙充滿水和油兩種流體，且均滿足零初始條件。本文算例選取的材料參數(shù)取自文獻[18－19]，具體數(shù)值如表1所示。

表1 材料參數(shù)Table 1 Material parameters

5.1 算例1

截斷項數(shù)N既關(guān)系到計算的精度，又直接影響計算的效率。為考慮N對計算結(jié)果的影響，選定上、下表面都完全透水、透氣，下表面固定，上表面作用單位階躍荷載作為邊界條件，即

圖2給出了邊界條件（74）時，兩種流體的動力滲透系數(shù)kw、kn分別為80、20 m4/(N·s)，N值分別為5、50、500、5 000時，土層中部ζ=0.5處總應(yīng)力響應(yīng)時程曲線。

由圖 2可以看到，N=5時，曲線不能完整反映土層的動力響應(yīng)。隨著N值的增加，計算響應(yīng)越來越精確。當 N= 5 000時，可以清晰的觀察到三類波的傳播，驗證了本文方法的合理性與有效性。綜合考慮計算效率和計算精度，如沒有特別說明，本文算例在階躍荷載作用時取 N= 5 000，正弦荷載作用時取 N= 1 000。

圖2 不同截斷項數(shù)下的響應(yīng)時程曲線Fig.2 Numerical results for different partial sums of the series solution

為進一步驗證本文方法，令非飽和多孔介質(zhì)內(nèi)的干相和濕相兩種流體都為水，此時，非飽和多孔介質(zhì)模型就退化為飽和多孔介質(zhì)。對于非飽和多孔介質(zhì)模型，邊界條件仍采用式（74）。對于飽和多孔介質(zhì)模型，采用如下邊界條件：

式中：p為飽和多孔介質(zhì)模型中的孔壓。

圖3、4給出了動力滲透系數(shù)為10-2m4/(N·s)時，分別用本文方法和文獻[10]方法計算得到ζ= 0.5處的孔壓和流體相對位移響應(yīng)結(jié)果。

由圖3、4可以看到，用兩種方法計算得到的孔壓結(jié)果幾乎完全一致，而飽和模型中的流體相對位移等于非飽和模型中濕相流體和干相流體的相對位移之和，從而進一步驗證了本文提出的非飽和多孔介質(zhì)動力響應(yīng)分析方法的有效性。

圖3 不同模型的孔隙流體壓力-時程曲線對比Fig.3 Pore pressure responses contrast of different models

圖4 不同模型的流體相對位移-時程曲線對比Fig.4 Fluid relative displacement responses contrast of different models

5.2 算例2

本算例同樣選取式（74）作為邊界條件，圖 5為不同動力滲透系數(shù)下土層ζ=0.5處濕相流體孔壓pw的響應(yīng)時程曲線。

由圖5可以看到，波在非飽和多孔介質(zhì)中的傳播與衰減特點與飽和多孔介質(zhì)中相似，但在非飽和多孔介質(zhì)中明顯存在 3類壓縮波。如圖 5(a)、5(b)所示，當動力滲透系數(shù)較大時，這3類壓縮波均清晰可見。大量計算表明，當流體滲透系數(shù)的數(shù)量級大于10-2時，多孔介質(zhì)的動力響應(yīng)幾乎完全一致；從 10-5量級開始，隨著滲透系數(shù)減小，流-固間的耦合作用逐漸增大，波的耗散逐漸加快，并在 10-8量級附近達到最大，此時3類壓縮波都迅速衰減，看不到明顯的波動現(xiàn)象。圖5(d)表明，當滲透系數(shù)很小時，固體和兩種流體之間存在很強的黏滯耦合作用，此時多相孔隙介質(zhì)呈現(xiàn)單相介質(zhì)的波動特性。

5.3 算例3

本算例采用如下邊界條件：

即土層上、下表面均完全透水、透氣，底部固定，土層頂部受到一個完整周期的正弦脈沖荷載作用。

圖6為不同動力滲透系數(shù)下土層ζ=0.2處的總應(yīng)力-時程曲線。從圖中可以更清晰、直觀地反映出不同滲透系數(shù)下非飽和土中壓縮波的傳播和衰減情況。在滲透系數(shù)較大時，流-固耦合作用很小，3類壓縮波的獨立傳播非常明顯，當滲透系數(shù)減小到10-8量級時，受流-固間的黏滯耦合作用，此時只能觀察到傳播速度較快的第1類波的波形，且波形不完整，說明此時波的耗散非?？臁．敐B透系數(shù)進一步減小到10-11量級時，只能觀測到1類波的波形，且波形是完整的，說明此時的多孔介質(zhì)性質(zhì)接近于單相介質(zhì)。

圖5 階躍荷載下ζ =0.5處濕相流體孔壓響應(yīng)曲線Fig.5 Pore pressure responses to step loading at ζ =0.5

圖6 正弦脈沖荷載作用ζ =0.2處總應(yīng)力響應(yīng)曲線Fig.6 Stress responses to sinusoidal impulse loading at ζ =0.2

5.4 算例4

本算例分析兩種流體的飽和度變化對非飽和多孔介質(zhì)動力響應(yīng)的影響，采用如下邊界條件：

動力滲透系數(shù)取kw=80 m4/(N·s)，kn=20 m4/(N·s)，這里無量綱過程中的m=λ+2μ+

圖7給出了 4種不同飽和度條件下，介質(zhì)層ζ= 0.8處的固體位移和兩種流體位移-時程曲線。

圖7 不同飽和度下ζ =0.8處位移響應(yīng)曲線Fig.7 Displacement responses for different saturations at ζ =0.8

由圖7可知，非飽和多孔介質(zhì)內(nèi)部兩種流體的不同比例會對介質(zhì)的瞬態(tài)響應(yīng)產(chǎn)生顯著影響。圖7(a)中可以清楚地觀察到 3類壓縮波先后傳到ζ=0.8處。飽和度的改變對傳播速度最快的第1類波速度影響很小，而對第2、2類波的波速影響較大。波速隨著濕相流體飽和度的增加而減小。此外，飽和度的改變對3類壓縮波的振動幅值和振動相位也有明顯影響。

6 結(jié) 語

本文基于 Zienkiewicz建立的非飽和多孔介質(zhì)波動方程，以一類典型的混合邊界條件為例提出了一種適用于任意邊界條件、任意材料參數(shù)的半解析方法。通過算例，驗證了本文所建立模型的正確性，探討了非飽和多孔介質(zhì)中波傳播的特性。

（1）非飽和多孔介質(zhì)內(nèi)存在三類壓縮波，波速最快的第一類波衰減較小，波速較慢的第二、三類波衰減較明顯。

（2）當兩種流體的動力滲透系數(shù)大于10-2m4/(N·s)量級時，介質(zhì)層的瞬態(tài)響應(yīng)基本上沒有差別；壓縮波的衰減隨著動力滲透系數(shù)的減小逐漸變快，并在10-8m4/(N·s)量級附近達到最大；當滲透系數(shù)小于10-10m4/(N·s)量級時，非飽和多孔介質(zhì)基本上呈現(xiàn)單相介質(zhì)的波動特性。

（3）三類壓縮波波速均隨著濕相流體飽和度的增加而減小。同時，當濕相流體飽和度的增加時，第一類壓縮波造成固體和流體的位移幅值均有所減少。第二類壓縮波造成的固體位移和濕相流體位移減少，而干相流體的位移增加。與之相反的，波速最慢的第三類壓縮波造成的固體位移和濕相流體位移增加，而干相流體的位移則有所減少。

[1]COUSSY O. Poromechanics[M]. West Sussex: John Wiley & Sons, 2004.

[2]BIOT M A. Theory of elastic waves in a fluid-saturated porous solid, I: Low-frequency range[J]. The Journal of the Acoustical Society of America, 1956a, 28(2): 168－178.

[3]BIOT M A. Theory of elastic waves in a fluid-saturated porous solid, II: High-frequency range[J]. The Journal of the Acoustical Society of America, 1956b, 28(2): 179－191.

[4]GARG S K, NAYFEH A H, GOOD A J. Compressional waves in fluid-saturated elastic porous media[J]. Journal of Applied Physics, 1974, 45(5): 1968－1974.

[5]SIMON B R, ZIENKIEWICZ O C, PAUL D K. An analytical solution for the transient response of saturated porous elastic solids[J]. International Journal for Numerical and Analytical Methods in Geomechanics,1984, 8(4): 381－398.

[6]GAJO A, MONGIOVI L. Analytical solution for the transient response of saturated linear elastic porous media[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1995, 19(6): 399－413.

[7]SCHANZ M, CHENG A H D. Transient wave propagation in a one-dimensional poroelastic column[J].Acta Mechanica, 2000, 145(1): 1－8.

[8]凌道盛, 張飛霞, 單振東, 等. 單層不可壓縮飽和多孔介質(zhì)一維瞬態(tài)響應(yīng)精確解[J]. 巖土力學(xué), 2011, 32(5):1033－1038.LING Dao-sheng, ZHANG Fei-xia, SHAN Zheng-dong,et al. Exact solution for one-dimensional transient response of single-layer incompressible fluid-saturated porous media under arbitrary vertical loadings[J]. Rock and Soil Mechanics, 2011, 32(5): 1033－1038.

[9]SHAN Z D, LING D S, DING H J. Exact solutions for one-dimensional transient response of fluid-saturated porous media[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2011, 35(4):461－479.

[10]凌道盛, 房志輝, 單振東, 等. 單層飽和多孔介質(zhì)一維瞬態(tài)響應(yīng)半解析解[J]. 巖石力學(xué)與工程學(xué)報, 2011,30(8): 1683－1689.LING Dao-sheng, FANG Zhi-hui, SHAN Zhen-dong,et al. A semi-analytical solution for one-dimensional transient response of single layered saturated porous media[J]. Chinese Journal of Rock Mechanics and Engineering, 1994, 34(2): 131－136.

[11]LI P, SCHANZ M. Wave propagation in a 1-D partially saturated poroelastic column[J]. Geophysical Journal International, 2011, 184(3): 1341－1353.

[12]SHAN Z D, LING D S, DING H J. Exact solutions for one-dimensional transient response of unsaturated single-layer porous media[J]. Computers and Geotechnics, 2013, 48: 127－133.

[13]BRUTSAERT W. The propagation of elastic waves in unconsolidated unsaturated granular mediums[J]. Journal of Geophysical Research, 1964, 69(2): 243－257.

[14]SCHANZ M, DIEBELS S. A comparative study of Biot′s theory and the linear theory of porous media for wave propagation problems[J]. Acta Mechanica, 2003, 161(3):213－235.

[15]ZIENKIEWICZ O C, CHAN A H C, PASTOR M, et al.Computational geomechanics with special reference to earthquake engineering[M]. Chichester: John Wiley &Sons, 1998.

[16]LI X, ZIENKIEWICZ O C, XIE Y M. A numerical model for immiscible two-phase fluid flow in a porous media and its time domain solution[J]. International Journal for Numerical Methods in Engineering, 1990, 30(6):1195－1212.

[17]鐘萬勰. 結(jié)構(gòu)動力方程的精細時程積分法[J]. 大連理工大學(xué)學(xué)報, 1994, 34(2): 131－136.ZHONG Wan-xie. On precise time-integration method for structural dynamics[J]. Journal of Dalian University of Technology, 1994, 34(2): 131－136.

[18]LO W C, SPOSITO G, MAJER E. Wave propagation through elastic porous media containing two immiscible fluids[J]. Water Resources Research, 2005, 41:W02025.

[19]ARNTSEN B, CARCIONE J M. Numerical simulation of the Biot slow wave in water-saturated Nivelsteiner sandstone[J]. Geophysics, 2001, 66(3): 890－896.