雙變量耦合作用對非飽和巖土波動特性的影響研究

胡亞元

( 浙江大學 濱海和城市巖土工程研究中心，杭州 310027)

巖土介質中波的傳播特性及動力響應研究一直受到理論力學界和工程界的重視。它不僅對彈性動力學發展具有重要理論價值，而且在地震學、地球物理學、石油工程、土木工程、動力機器基礎和防護工程等工程領域具有廣闊應用前景。近年來，隨著非飽和巖土力學的發展，對波在非飽和巖土的傳播特性和動力響應的研究越來越深入，但采用雙變量理論(雙應力或雙應變)研究非飽和巖土波動特性的文獻目前還極其少見，因此有必要開展雙變量理論框架下的非飽和巖土波動特性研究。

建立非飽和巖土波動理論控制方程的方法主要有兩種。一種是Biot創建的工程法[1-2]。Brutsaert[3]和Santos等[4]把非飽和巖土看成由固液氣三相連續組成的多孔介質，基于Lagrangian坐標系把Biot飽和巖土波動理論擴展到非飽和巖土波動理論，理論證明在非飽和巖土中除與飽和巖土對應的剪切波S波、快縱波P1波和第一慢縱波P2波外，還存在與第二種流體對應的第二慢縱波P3波。 蔡袁強等[5]和徐明江[6]從單變量理論出發，應用Skempton有效應力公式[7]，建立了非飽和巖土波動方程并深入研究非飽和巖土四種體波(包括三個縱波和一個橫波)的傳播特性。另一種是混合物理論法。Vardoulakis等[8]采用混合物理論建立了非飽和巖土波動方程，徐長節等[9]采用Vardoulakis波動方程深入研究了非飽和土中三種彈性縱波的波動特性。Wei等[10-11]根據相分離原理建立了非飽和巖土波動控制方程，黃義和張引科[12]在忽略固相和液相材料變形的條件下根據混合物理論建立了非飽和土的波動控制方程，系統建立了求解非飽和土動力響應的邊界元理論。Chen等[13]根據Wei等建立的非飽和巖土波動方程，忽略流體組分體積分數對流體自由能影響，深入分析了非飽和巖土的表面波傳播特性。上述這些研究深化了工程界對彈性波在非飽和巖土介質中傳播特性的認識，極大地推動了非飽和巖土波動理論的發展。但由于非飽和土中固液氣三相相互作用的復雜性，對非飽和巖土波動控制方程的研究仍然存在一些不足，主要表現在兩個方面:①某些非飽和巖土波動方程中的模型參數取值比較困難，比如徐明江指出，按照混合物理論建立的一些非飽和巖土波動方程存在這類困難。②Fredlund理論[14]表明，非飽和巖土的力學特性需要采用雙變量才能得到合理地反映,因此非飽和巖土的波動方程也需要依據雙變量理論來建立。然而目前采用雙變量理論來研究非飽和巖土波動特性和動力響應的文獻還很少見，亟需填補這方面的研究空缺。

本文從Houlsby等[15]，Borja[16]和胡亞元[17]創建的非飽和多孔介質工程力學理論出發，根據能量平衡方程和熱力學局部平衡原理構建自由能勢函數，在假定各組分材料之間以及它們與骨架之間的力學特性均不相互耦合的前提下建立了以雙變量理論為基礎的非飽和巖土彈性本構方程。根據不同雙變量建立的線彈性本構方程相互一致的性質，推導出彈性本構方程模型參數與非飽和巖土常規土工參數之間的換算公式。利用非平衡態熱力學熵產公式構建滲流耗散率勢函數，獲得液氣兩相的達西滲流本構方程。把上述本構方程與動量守恒方程相結合，獲得一個以雙應力變量為基礎的非飽和巖土線彈性波動控制方程。最后，通過算例分析了雙應力變量耦合效應對非飽和土三個縱波和一個橫波波速的影響規律。

1 波動方程

1.1 動量和動量矩守恒

Sr=nL/(nL+nG)

(1)

根據各組分體積分數的定義，有：

nS+nL+nG=nS+n=1

(2)

(3)

設σα為第α組分的柯西應力，σ為混合物總柯西應力，根據混合物理論有：

(4)

(5a)

(5b)

不考慮力偶供應量，根據各相的角動量守恒方程可得各相的柯西應力張量是對稱張量。

1.2 能量平衡方程

(6)

定義Bishop 平均有效應力為[15-16，19-20]：

(7)

(8)

(9)

令Pr=-σ∶1/3為多孔介質的總壓力，則定義固相的真實壓力為：

(10)

(11)

(12)

(13)

設εS為ES在小應變條件下的固相應變張量，則根據ES和EH的關系有：

(14a)

根據組分體應變與其密度之間存在的數學關系[21],由ρL=SrnρRL，ρG=(1-Sr)nρRG并利用式(12)可得小應變條件下液氣兩相的組分體應變εLv和εGv分別為:

(14b)

(14c)

(15)

根據大應變和小應變之間物理量的關系，由式(11)～(13)和式(15)可得：

(16)

式(16)即是小應變條件下非飽和多孔介質混合物的能量平衡方程。

1.3 線彈性方程和線性滲流方程

根據能量平衡方程，等存性公設和Lagrangian乘子法理論，彈性條件下非飽和多孔介質的自由能可假定為：

E(η,εH,Sr,εRα,λL,λG)=

(17)

對式(17)進行微分得：

(18)

根據熱力學局部平衡假定，可得：

(19a)

(19b)

引入Helmhotlz自由能Ψ*(θ,εH，Sr，εRα)，它與內能之間的關系為：

Ψ*(θ,εH，Sr,εRα)=E*(η,εH,Sr,εRα)-θη

(20)

把式(20)兩邊微分再把式(19a)代入得：

(21)

(22)

式中：DHH、KRαγ、Krr、Kθθ、KRHα、KHr、KHθ、KRrα、KRθα、Kθr是模型彈性系數，KRαγ=KRγα，α∈{S,L,G}，γ∈{S,L,G}。把式(22)代入到式(21)得：

(23a)

(23b)

(23c)

(23d)

在目前通行的飽和和非飽和巖土力學中，均假定各組分材料的本構模型在形式上與其單相時的本構模型相同[14,20]，即其壓力只與其本身的體積變形和溫度有關而與其它相的變形無關。在該假定下，各組分的材料變形之間以及它們與骨架變形之間均相互解耦。同時假定波動過程處于熱力學等溫過程，故取θΔ=0。利用這些假定，式(23a)～(23d)簡化為：

(24a)

(24b)

nαPα=KRααεRαα∈{S,L,G}

(24c)

式(24a)～(24b)表明，當Bishop平均應力需要采用雙應變變量表達時，修正吸力也要采用雙應變變量來表示。如果修正吸力不相應地采用雙應變變量來表示，則由此獲得的波動方程剛度矩陣就會失去對稱性。同時，由于土工試驗表明巖土和混凝土等非飽和多孔介質的基本力學性質及其本構關系需要采用雙變量(雙應變或雙應力)來表示，因此式(24a)～(24c)是其最簡單形式。從式(20)～(24c)的推導過程看，與文獻[12]相比，本文不但考慮了固相和液相自身的材料變形，而且還大大減少了模型參數，簡化了獲得非飽和巖土線彈性本構關系的推導過程。與蔡袁強等基于工程法提出的非飽和巖土單變量線彈性本構模型相比，本文式(24a)～(24b)考慮了骨架變形和飽和度之間的耦合作用，它不僅與當前非飽和巖土力學的雙變量理論之間更加一致，而且因為修正吸力也采用了雙應變變量來表示，所以滿足彈性功互易定理，線彈性本構方程的剛度矩陣相應地也成為對稱矩陣。

當固相骨架為各向同性材料時，式(24a)～(24b)簡化為：

(25a)

(25b)

式中：λH和G為骨架Lame常數；KHr為骨架和毛細吸附作用之間的耦合剛度；Krr為修正吸力關于飽和度的剛度；1為二階單位張量；14為四階單位張量，(14)ijkl=(δikδjl+δilδjk)/2。對式(25a)～(25b)以及式(24c)求逆得：

(26a)

(26b)

(26c)

式中：EH為骨架的楊氏模量；υ有骨架的泊松比；EHr為骨架和毛細吸附耦合作用之間的彈性模量；Err為飽和度隨吸力變化的彈性模量。由于式(25a)～(25b)和式(26a)～(26b)互逆，固有：

(27a)

(27b)

(27c)

(28a)

(28b)

(28c)

把式(28a)～(28c)代入到式(26c)得：

(29a)

(29b)

把式(26a)，式(28a)～(28c)和式(29a)代入式(14a)得：

(30)

(31a)

(31b)

把式(26a)～(26c)、式(28a)～(28c)和式(29a)～(29b)代入式(31a) ～(31b)得：

(32a)

(32b)

由于在波動理論中不考慮熱流量和外熱供應，故qα=0，rα=0。把式(18)代入到式(16)并利用式(19a)～(19b)得：

(33)

(34)

為了獲得液氣兩相的線性化達西定律，把式(34)中的耗散率勢函數D*表達為:

(35)

根據式(34)有(β∈{L,G})：

(36)

從上面的推導過程可以看出，由于本文基于多孔介質工程力學來建立非飽和巖土的線彈性方程，故方程中的各個參數均有明確的物理含義，這就為建立模型參數和常規土工參數之間的換算公式提供了方便。

1.4 模型參數與土工參數之間的關系研究

在巖土力學中，一般都假定組分材料的本構關系與單相時的本構關系相同。令第α組分材料單相時與壓強Pα對應的體積模量為Kα，則有：

KRαα=nαKα

(37)

在非飽和巖土土工試驗中，保持PL和PG恒定，變化σ，即完全排水和排氣，此時可以測定多孔固體的體積模量Kb和剪切剛度G，由式(30)和式(37)可得：

(38a)

(38b)

同時，保持總應力σ和PG恒定，變化PL，可以測定多孔固體體積應變隨吸力變化的彈性模量ESr和土-水特征曲線方程Sr=Sr(s)，根據式(30)和式(38a)有：

(39a)

根據總應力σ和氣體孔壓PG恒定下的土-水特征曲線方程Sr=Sr(s)、式(26b)和式(39a)有：

(39b)

令：

(40a)

(40b)

式中：φ為Biot系數，χ為Bishop系數。把式(38a)～(38b)、式(39a)和式(40a)～(40b)代入到式(30)和式(32a)～(32b)并利用εSv=εS∶1和σm=σ∶(1/3)得：

(41a)

(41b)

(41c)

令：

(42a)

(42b)

(42c)

cSL=φ[χcLL+(1-χ)cLG]

(42d)

cSG=φ[χcLG+(1-χ)cGG]

(42e)

(43a)

M=φχcSL+φ(1-χ)cSG=

φ2[χ2cLL+2χ(1-χ)cLG+(1-χ)2cGG]

(43b)

(43c)

求解式(41a)～(41c)并利用式(42a)～(42e)和式(43a)～(43b)得：

(44a)

(44b)

根據σ=(σ-σm1)和式(44a)并利用式(43c)得：

(45a)

(45b)

2 波速公式

(46a)

(46b)

(47a)

(47b)

設液氣兩相的滲透系數分別為kLL和kGG，液氣兩相的耦合滲透系數為kLG和kGL，則根據bLL、bLG和bGG與kLL、kLG、kGG和之間的互逆關系可得：

(48a)

(48b)

式(48a)中ηL和ηG為液氣兩相的黏滯系數。

令eS=εS∶1，把式(45a)～(45b)代入到式(47a)～(47b)得：

(49a)

(49b)

由于體力只與波動方程的特解有關，其本身不影響體波波速特性，故為簡便起見略去式(49a)～(49b)中的體力項，對式(49a)～(49b)求散度并利用ξβ=-▽·wβ得：

(50a)

(50b)

顯然，式(50a)～(50b)反映的是壓縮波的傳播特性，令

{eS,ξL,ξG}={AS,AL,AG}exp[j(ωt-lp·r]

(51)

(52)

式中：

(53a)

(53b)

要使式(52)有非零解，其系數行列式必須等于零，故有：

(54)

(55)

理論分析有時把式(54)展開較為方便，由此可得：

(56)

式中:

(57a)

(57b)

(57c)

(57d)

求解式(56)同樣可獲得非飽和土的三個縱波波數。

略去式(49a)～(49b)中的體力項，對式(49a)～(49b)求旋度得：

(58a)

(58b)

式(58a)～(58b)反映的是剪切波的傳播特性，令

▽×{uS,wF,wG}=

{Bs,BL,BG}exp[j(ωt-lS·r)]

(59)

式中：ω為頻率，r為位移矢,lS為波矢量，lS為波數(波矢量模)。把式(59)代入到式(58a)～(58b)得：

(60)

式(60)中：

(61)

故剪切波的相速度和特征衰減系數為：

(62)

3 算例分析雙應力耦合效應對波速特性的影響

本節采用數值方法分析雙應力耦合效應對非飽和土四個體波波速的影響。模量ESr反映了多孔固體和毛細吸附之間的耦合效應，它與χ之間的關系見式(40a)。為了反映非飽和土雙應力變量之間的耦合效應程度，定義雙應力變量的耦合效應系數為：

(63)

本節數值分析時土-水特征曲線采用van Genuchten模型，其表達式為[22]：

Sr=[1+(αvgS)nvg]-mvg

(64)

式中：αvg、mvg和nvg是van Genuchten模型的三個參數。相應的滲透系數kLL、kLG、kGL和kGG為：

(65a)

(65b)

kGL=kLG=0

(65c)

式中：κ為本征滲透系數。根據式(39b)、式(40a)和式(64)有：

(66)

根據式(48a)-(48b)和式(65a)～(65c)有：

(67)

本文選用非飽和土作為算例的研究對象，故固體是土體。數值分析所用的固、液和氣的材料彈性模量按常規取值；固、液和氣的物理參數、土體骨架的彈性模型參數和van Genuchten模型參數按文獻[20]提供的土工參數換算獲得，如表1所示。

表1 非飽和土的物理力學參數Tab.1 Physical and mechanical parameters of unsaturated soil

上述規律的力學機理在于：在孔隙流體中，非飽和土的P1和P2縱波波速在一般情況下較大地受KL/nL和KG/nG這兩個數值中的較小者影響。由于液體的體積模量KL遠大于氣體模量KG，當飽和度小于90%時，KG/nG比KL/nL小的多，P1和P2縱波波速主要受氣體力學參數控制。當飽和度大于90%時，隨著KL/nL與KG/nG數值逐漸接近，液體KL/nL對P1和P2縱波的影響越來越明顯，當飽和度接近100%時，KG/nG趨向于無窮大，KL/nL成為較小值，此時P1和P2縱波波速主要受液體力學參數控制。正如學者徐明江和蔡袁強等根據數值分析成果所指出的那樣，當飽和度大于90%時，P1和P2縱波從主要受空氣體積模量控制過渡到主要受液體體積模量控制，因此，圖1中P1和P2縱波波速在這個飽和度區間的變化幅值比較大。固流兩相的滲透黏滯效應對P2縱波具有制約作用，故P1縱波波速的變化幅值比P2縱波更加明顯。對于P3縱波，當飽和度等于0%時，孔隙中只有氣體一種流體，它屬于氣體飽和土；當飽和度等于100%時，孔隙中變為只有液體一種流體，它屬于液體飽和土，由于它們都是飽和土，這兩種極端情況下的P3縱波波速均等于零。因此當飽和度從零變化到100%時，P3縱波波速的變化規律是從零開始隨著飽和度逐漸增加，到極值后再逐漸減小回歸到零。與P1和P2縱波一樣，這兩個階段的轉折點也在飽和度90%左右。P3縱波波速比較小，只有每秒零點幾米。由于基數比較小，即使增加量比較小，相對影響值也比較大。這是當飽和度小于90%時圖1中P3波波速變化影響較為明顯的原因。當飽和度大于90%時，由于非飽和土逐漸逼近于液體飽和土，第三縱波的波速重新回歸為零，這是P3波當飽和度大于90%時波速變化規律較為明顯的原因。

(a) P1 波

(b)P2 波

(c) P3 波

(d) S 波圖1 非飽和土中四種體波隨飽和度變化圖Fig.1 the velocities of four body waves vary with the saturated degree in unsaturated soil

圖1表明當飽和度大于95%時，三種縱波的波速變化劇烈。特別是P1波，當飽和度大于95%后，波速遠大于飽和度小于90%時的波速。本文重點是研究雙應力變量耦合效應對非飽和土四種體波波速的影響，當飽和度大于95%后非飽和土四種體波波速的劇烈變化會干擾雙應力變量耦合效應對非飽和土四種體波波速影響的分析，因此圖2和圖3中飽和度的取值范圍調整為(0.05～0.95)。雖然沒有給出飽和度大于95%后雙應力變量耦合效應影響體波波速的關系曲線，但根據圖2和圖3中體波的發展趨勢可以容易預測其影響規律。

圖2給出了頻率20 Hz時雙應力變量耦合效應系數對四種體波的影響，從圖2中可以看出，雙應力變量耦合效應對剪切波波速沒有影響，這也可以從公式(60)～(61)中看出。雙應力變量耦合效應對P2波和P3波的影響很小，可以忽略不計。當飽和度較小時，雙應力變量耦合效應對P1波波速有明顯的影響，雙應力變量耦合效應系數越大，對P1波波速的影響越明顯。但影響程度隨著飽和度的增大迅速減小。上述分析表明，當飽和度較小時，雙應力變量耦合效應對P1波速的影響是不可忽略的。

(a) P1 波

(b)P2 波

(c) P3 波

(d) S 波圖2 頻率20Hz時四種體波波速隨系數變化圖Fig.2 the velocities of four body waves vary with the coefficient at the frequency 20 Hz

圖3呈現了頻率200 Hz時雙應力變量耦合效應系數對四種體波的影響規律，所反映的規律特征與圖2一致。雙應力變量耦合效應對剪切波波速沒有影響，對P2波和P3波的影響很小，可以忽略不計。當飽和度較小時，雙應力變量耦合效應對P1波波速有明顯的影響，但影響程度隨著飽和度的增大迅速減小。通過對比圖2和圖3還可以發現，在低頻段，頻率對P1波波速和S波波速影響較小，但對P2波和P3波的影響較大。但頻率從20 Hz增加到200 Hz時，P2波波速增加了約3倍，P3波的波速增加了約3.2倍。因此，即使在低頻段，頻率對P2波和P3波波速的影響也無法忽略。

(a) P1 波

(b)P2 波

(c) P3 波

(d) S 波圖3 頻率200 Hz時四種體波波速隨系數變化圖Fig.3 the velocities of four body waves vary with the coefficient at the frequency 200 Hz

4 結 論

本文根據非飽和多孔介質工程力學理論，建立了基于雙變量理論的非飽和巖土波動理論，獲得了以下研究成果：

(1)本文所依據的非飽和多孔介質工程力學理論以巖土力學中常用的應力應變物理量作為建模所用的狀態變量，由它推導出的波動方程不但滿足彈性力學的互易定理，從而使剛度矩陣具有對稱性要求，彌補了當前根據工程法獲得的非飽和巖土波動方程剛度矩陣不對稱的缺陷；而且容易根據常規土工試驗選取模型參數，從而更有利于工程應用。

(2)與不考慮雙變量耦合效應的波速傳播特性相比，雙應力變量耦合效應對非飽和巖土的剪切波波速沒有影響。對P2波和P3波的影響很小，可以忽略不計。對P1波波速的影響在飽和度較小時比較明顯，耦合效應越大，對P1波波速的影響越明顯。但其影響程度隨著飽和度的增大迅速減小。

(3)通過對比圖2和圖3還可以發現，在低頻段，頻率對P1波波速和S波波速影響較小，但對P2波和P3波的影響較大。但頻率從20 Hz增加到200 Hz時，P2波波速增加了約3倍，P3波的波速增加了約3.2倍。因此，即使在低頻段，頻率對P2波和P3波波速的影響也無法忽略。

參 考 文 獻

[1] BIOT M A.Theory of propagation of elastic waves in a fluid-saturated porous solid (I.low-frequency range) [J].J.Acoust.Soc.Am., 1956, 28(2): 168-178

[2] BIOT M A.Theory of propagation of elastic waves in a fluid-saturated porous solid(Ⅱ.higher-frequency range) [J].J.Acoust.Soc.Am., 1956, 28(2): 179-191

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

[4] SANTOS J E, RAVAZZOLI C L, GAUZELLINO P M, et al.Simulation of waves in poro-viscoelastic rocks saturated by immiscible fluids-Numerical evidence of a second slow wave[J].J.Comput.Acoust., 2004, 12(1): 1-21.

[5] 蔡袁強, 李保忠, 徐長節.兩種不混溶流體飽和巖石中彈性波的傳播[J].巖石力學與工程學報, 2006, 25(10): 2009-2016.

CAI Yuanqiang, LI Baozhong, XU Changjie.Analysis of elastic wave propagation in sandstone saturated by two immiscible fluids [J].Chinese Journal of Rock Mechanics and Engineering, 2006, 25(10): 2009-2016.

[6] 徐明江.非飽和土地基與基礎的動力響應研究[D].廣州：華南理工大學, 2010.

[7] SKEMPTON A W.Effective stress in soils, concrete and rocks[J].Butterwoths: London UK, 1961: 4-16.

[8] VARDOULAKIS I, BESKOS D E.Dynamic behavior of nearly saturated porous media[J].Mechanics of Materials, 1986(5): 87-108.

[9] 徐長節, 史焱永.非飽和土中波的傳播特性[J].巖土力學, 2004, 25(3): 354-358.

XU Changjie, SHI Yanyong.Characteristics of wave propagation in unsaturated soils[J].Rock and Soil Mechanics, 2004, 25(3): 354-358.

[10] WEI C F, MURALEETHARAN K K.A continuum theory of porous media saturated by multiple immiscible fluid: I.Linear poroelasticity[J].International Journal of Engineering Science, 2002, 40:1807-1833.

[11] WEI C F, MURALEETHARAN K K.A continuum theory of porous media saturated by multiple immiscible fluids: II.Lagrangian description and variational structure[J].II.International Journal of Engineering Science,2002, 40: 1835-1854.

[12] 黃義，張引科.非飽和土本構關系的混合物理論[J].應用數學和力學, 2003, 24(2): 111-137.

HUANG Yi, ZHANG Yinke.Mixture theory of constitutive relation of unsaturated soil[J].Applied Mathematics and Mechanics, 2003, 24(2): 111-137.

[13] CHEN W Y, XIA T D, HU W T.A mixture theory analysis for the surface-wave propagation in an unsaturated porous medium[J].International Journal of Solids and Structures, 2011, 48(16): 2402-2412.

[14] 弗雷德隆德 D G, 拉哈爾佐 H.非飽和土土力學[M].陳仲頤, 張在明, 陳愈炯, 譯.北京: 中國建筑工業出版社, 1997.

[15] HOULSBY G T, PUZRIN A M.Principles of Hyper- plasticity[M].Springer: London, 2006:211-239.

[16] BORJA R I.On the mechanical energy and effective stress in saturated and unsaturated porous continua[J].International Journal of Solids and Structures, 2006, 43: 1764-1786.

[17] 胡亞元．非飽和多孔巖石的熱力學本構理論研究[J].浙江大學學報工學版, 2017,51(2): 255-263.

HU Yayuan, Thermodynamics-based constitutive theory for unsaturated porous rock[J].Journal of Zhejiang University (Engineering Edition), 2017,51(2): 255-263.

[18] BOWEN R M.Theory of mixture.In AC Eringen (ed): Continuum Physics, Vol III[M].Academic Press: New York, San Francisco, London, 1976

[19] BISHOP A W.The principle of effective stress[J].Teknisk Ukeblad, 1959, 106(39): 113-143.

[20] WHEELER S J, SHARMA R S, BUISSON M S R.Coupling of hydraulic hysteresis and stress-strain behavior in unsaturated soil[J].Geotechnique, 2003,53(1): 41-54.

[21] 黃筑平.連續介質力學基礎[M].北京:高等教育出版社, 2003:78-122.

[22] VAN GENUCHTEN M T.A closed form equation for predicting the hydraulic conductivity of unsaturated soils[J].Soil Science Society of America Journal, 1980, 44(5):892-898