非飽和雙重孔隙介質的勢函數本構方程及應用

胡亞元

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

隨著對巖土孔隙尺度分布特征的深入了解，發現許多巖土材料具有明顯的雙重孔隙結構特征[1-2]。巖土工程界把較小尺寸的孔隙仍稱為孔隙，把較大尺寸的孔隙稱為裂隙，把具有雙重孔隙結構的巖土材料統稱為雙重孔隙介質。混合物理論是非飽和雙重孔隙介質常用的建模方法之一[3-8]。Borja等[3]推導了非飽和雙重孔隙介質的能量守恒方程，并根據功共軛量的力學性質提出采用Skempton型有效應力、孔隙內修正吸力和孔隙間平均修正吸力作為應力狀態變量來構建非飽和雙重孔隙介質的本構模型。Zhang等[4]基于Borja等[3]的建議，創建了非飽和雙重孔隙土的橫觀各向同性線彈性方程。Li等[5-6]假定固相和吸附水之間存在質量交換，通過變形功公式中的功共軛對來選擇應力和應變變量，建立非飽和雙重孔隙膨脹土的彈塑性模型，較好地揭示了Alonso等[1]試驗呈現出的變形和持水特性。Sanchez等[7]把土應變拆分成兩個獨立的固相應變之和，通過組合這兩個固相應變的本構模型來建立非飽和雙重孔隙土的彈塑性模型。Guo等[8]把Borja等[3]建立的能量守恒方程轉化為采用兩個固相應變和兩個有效應力相共軛表示的形式，建立了非飽和雙重孔隙土的雙有效應力彈塑性模型。上述研究有力地推動了非飽和雙重孔隙介質力學理論和本構模型的發展。

在非飽和多孔介質中，固液氣相的賦存空間均與孔隙狀態密切相關。孔隙變形同時影響固液氣三相的應變量，協調、傳遞和制約流固之間的耦合作用，在多相耦合機制中處于舉足輕重的關鍵紐帶地位[9]。然而在以往的混合物理論研究中，卻沒有把孔隙變形當做一個獨立的應變量來對待，從而忽略了孔隙變形與各相應變之間的內在關系，難以揭示多孔介質流固耦合作用的力學本質，限制了多孔介質力學理論的進一步深入發展。胡亞元[9]把孔隙變形作為一個獨立的應變量來對待，較好地解決了非飽和單重孔隙介質流固兩相本構關系之間的耦合問題。非飽和雙重孔隙介質存在兩種尺度的孔隙類型，需要兩個反映孔隙變形的應變來度量。為了實現這一目的，本文把雙重孔隙介質視為兩個單重孔隙介質的嵌套疊加。如在雙重孔隙土中，內含孔隙的團聚體是一個單重孔隙介質；把內含孔隙的團聚體整體視為基質，把團聚體間裂隙視為單重孔隙，又構成單重裂隙介質。雙重孔隙土可視為單重裂隙介質的基質中嵌套單重孔隙介質而成。又例如，孔隙-裂隙巖體可視為由含有單重孔隙的完整巖塊作為基質，以裂隙作為單重裂隙的裂隙巖體。前一個單重孔隙介質(含孔隙的完整巖塊)作為基質嵌套在后一個單重孔隙介質(裂隙巖體)之中。根據上述嵌套思路，本文首先把固相變形拆分為固相基質變形和每個孔隙率變化引起的固相骨架應變之和，然后從混合物理論出發推導非飽和雙重孔隙介質的能量平衡方程，建立了非飽和雙重孔隙介質的一般自由能勢函數本構理論框架。

線彈性模型是巖土工程固結和波動理論中必不可少的本構模型，如陳正漢[10]、Fredlund等[11]、周萬歡等[12]、胡亞元[13]和Moradi等[14]運用線彈性本構模型建立了非飽和單重孔隙土的波動和固結理論，用于建筑場地的波動和固結特性分析。Berryman等[15]、Khalili等[16]和Yang等[17]利用線彈性本構模型研究了飽和雙重孔隙場地的波動和固結特性。本文的另一項研究工作是從自由能勢函數本構方程出發，推導了非飽和雙重孔隙介質的線彈性本構方程，并運用它建立了膨潤土的固結控制方程。根據本文研究成果，能夠較方便地建立非飽和雙重孔隙介質波動和固結方程。

1 質量、動量和動量矩守恒方程

1.1 雙重孔隙介質的組分占比指標

雙重孔隙介質由固相基質(在土力學中也稱為土顆粒)，孔隙(或小孔隙)和裂隙(或大孔隙)組成，在孔隙和裂隙中各含有液體和氣體。設固相用S表示，孔隙用P表示，裂隙用F表示，液相用L表示，氣相用G表示。PL和PG分別表示孔隙中的液相和氣相；而FL和FG分別表示裂隙中的液相和氣相。把裂隙視為廣義的孔隙，令β為孔隙類別指標，β∈{P,F}，當β=P時為孔隙，而β=F時為裂隙；γ為流相類別指標，γ∈{L,G}，當γ=L時為液相，而γ=G時為氣相。φ表示體積分數，有：

(1)

1.2 質量守恒

根據混合物理論，把混合物各組分按體積分數平均到整個混合物空間后，固相、裂隙和孔隙中液相和氣相共同連續地占有非飽和雙重介質混合物的空間位置。令α為組成非飽和雙重孔隙介質的組分指標，α∈{S,PL,PG,FL,FG}，設第α組分的初始位置為Xa，在當前t時刻的空間位置為x，則每一組分的運動方程可表示為x=xα(Xα,t)，速度可表示為vα=dxα(Xα,t)/dt，加速度可表示為aα=d2xα(Xα,t)/dt2。對于定義在x和t上的標量場或矢量場Γα，基于α組分的物質導數的定義為

(dαΓα/dt)=(?Γα/?t)+gradΓα·vα

(2)

假定固相、液相和氣相之間不存在質量交換，但裂隙液相與孔隙液相以及裂隙氣相與孔隙氣相間存在質量交換，各組分的質量守恒方程為：

(dSρS/dt)+ρSdivvS=0

(3)

(dβγρβγ/dt)+ρβγdivvβγ=cβγ

(4)

式中β∈{P,F}和γ∈{L,G}，下同；cβγ表示β孔隙中γ流相的質量交換率，有

(5)

把固相作為非飽和雙重孔隙介質混合物的參考構形，令β孔隙中γ流相相對固相的擴散速度為Wβγ=vβγ-vS。把它和ρβγ=φβγρRβγ代入式(4)得

Wβγgradφβγ-(cβγ/ρRβγ)=0

(6)

(7)

(8)

(9)

把式(8)減去式(9)得

(10)

1.3 動量和動量矩守恒

令σ為非飽和雙重孔隙介質混合物的Cauchy總應力張量，σS為固相組分的Cauchy應力張量，σβγ為β孔隙中γ流相組分的Cauchy應力張量，有

(11)

令I為二階單位張量，注意到在混合物中應力以拉為正而壓力以壓為正，故非飽和雙重孔隙介質混合物上的總壓力為PT=-σ∶I/3；固相基質所受的真實壓力為PS=-σS∶I/(3φS)，β孔隙中γ流相的孔壓為PβγI=-σβγ/φβγ，由式(11)得

(12)

圖1給出了非飽和雙重孔隙介質中孔隙和裂隙介質的結構層次關系和應力關系。圖1(a)為非飽和雙重孔隙介質單元體，圖1(b)顯示了非飽和雙重孔隙介質總壓力與各組分壓力關系式(12)的力學內涵。

圖1 非飽和雙重孔隙介質特征單元體示意Fig.1 Schematic diagram of representative volume element of unsaturated double-porosity media

把作為裂隙介質基質的孔隙介質視為一個單獨的混合物來進行分析。非飽和孔隙介質的應力等于組成非飽和孔隙介質的各組分應力之和：

(13)

(14)

圖1(c)顯示了從非飽和雙重孔隙介質隔離出來的非飽和孔隙介質單元體，圖1(d)圖顯示了非飽和孔隙介質總壓力和各組分壓力關系式(14)的力學內涵。

(15)

(16)

假定第α組分的動量矩供應量為0，則由動量矩平衡方程可得σα應力張量是對稱張量。

2 能量平衡方程

(17)

(18)

(19)

(20)

(21)

式(20)～(21)表明，DC與孔隙介質在裂隙介質中的體積分數有關，孔隙介質是裂隙介質的基質，在裂隙介質中起到骨架作用，故DC稱為裂隙骨架變形率。同理DH稱為孔隙骨架變形率。知道DC和DH的力學內涵后，式(19)所表示的能量守恒方程的物理內涵就變得十分明顯。式(19)等式右邊的前兩項表示的是與裂隙骨架及吸力相關的變形能，第三四項(即第一個括號內的那兩項)表示的與孔隙骨架及吸力相關的變形能，第五項為各組分基質引起的變形能，第六七項是非飽和雙重孔隙介質各組分內部相互作用引起的能量耗散，最后三項為熱傳遞和熱交換引起的能量變化。

從式(19)推導過程可知，固相變形率DS被分解成三部分：DC、DH和固相基質體應變率dS?S/dt。在單重裂隙介質中，裂隙介質基質的體積分數和裂隙體積分數之和等于1，故DC也與裂隙介質中裂隙體積分數的變化相關。同理DH也與孔隙介質的孔隙體積分數變化相關。這種變形分解方式有利于突出孔隙變形在多孔介質力學多場耦合機制中的紐帶作用。固相應變的上述嵌套分解結果也可以采用連續介質力學變形梯度的極分解理論來解釋。如圖2所示，把雙重孔隙介質視為以孔隙介質為基質的裂隙介質，故固相變形梯度可分解為孔隙介質變形梯度和裂隙骨架變形梯度的乘積；再把孔隙介質視為由固相基質與孔隙組成，故孔隙介質變形梯度可分解為固相基質變形梯度與孔隙骨架變形梯度的乘積。

圖2 雙重孔隙介質變形梯度示意Fig.2 Schematic diagram of deformation gradient of double-porosity media

(22)

3 小應變線彈性本構方程

在小應變條件下，令εS為固相應變張量ES的近似值，εC為裂隙骨架應變張量UC的近似值，εH為孔隙骨架應變張量EH的近似值，有：

εS=εC+εH-(?S/3)I

(23)

εC=[(?uC/?XS)+(?uC/?XS)T]/2

(24)

εH=[(?uH/?XS)+(?uH/?XS)T]/2

(25)

?S=ln(ρRS/ρRS0)≈(ρRS-ρRS0)/ρRS0

(26)

εVS=εVC+εVH-?S

(27)

注意到dεVC/dt=DC∶I， dεVH/dt=DH∶I。根據式(20)～(21)并在小應變條件下略去高次項得：

εVC=-ln(φSP/φSP0)≈(φSP0-φSP)/φSP0

(28)

(29)

?βγ=ln(ρRβγ/ρRβγ0)≈(ρRβγ-ρRβγ0)/ρRβγ0

(30)

(31)

(32)

(33)

(34)

式中ρβγ0和φβ0分別是β孔隙中γ流相的初始密度和初始體積分數，φSO是固相的初始體積分數，SrPO和SrFO分別是孔隙和裂隙的初始飽和度。從式(27)、式(31)～(34)可以看出，εVS和εVFL、εVFG的計算式中均含有裂隙骨架體應變εVC；εVS和εVPL、εVPG的計算式中均含有裂隙和孔隙骨架體應變εVC和εVH，故εVS與εVFL、εVFG、εVPL、εVPG之間存在耦合作用，它們是通過εVC和εVH傳遞的。本文選用裂隙和孔隙變形引起的εVC和εVH作為獨立變量，可以把非飽和雙重孔隙介質的耦合機理通過顯式表達式呈現出來，這是以往的研究中從來沒有揭示過的。定義β孔隙中γ流相的滲入量為ζβγ=φβγ0(εVS-εVβγ)，利用式(27)和式(31)～(34)可得ζβγ的表達式為：

(35)

(36)

(37)

(38)

(39)

能量平衡方程根據小應變條件和式(22)近似等于

(40)

假定受力過程滿足熱力學局部平衡假定，則內能可表示為ξ=ξ(η,εC,εH,Srβ,?α)，對ξ求全微分后由式(40)和理性力學的Coleman關系可得：

(41)

(42)

(43)

引入Helmhotlz自由能ψ=ψ(θ,εC,εH,Srβ,?α)，它與內能之間的關系為ψ=ξ-θη，對ψ求全微分后由式(41)～(42)得：

(44)

(45)

式(44)～(45)中的ψ是自由能勢函數，就熱力學性質而言是可逆的，故式(44)～(45)是非飽和雙重孔隙介質的一般彈性本構關系。非飽和雙重孔隙介質的粘塑性力學特性可以根據非平衡態熱力學從式(43)獲得，相關研究可借鑒文獻[9]，受篇幅限制本文不擬展開研究。

(46)

式中KCC、KHH、Krβζ、Kαχ、KCH、KCrβ、KCα、KHrβ、KHα和KJβα為模型的彈性參數，Kαχ=Kχα，Krβζ=Krζβ，α,χ∈{S,PL,PG,FL,FG}，β,ζ∈{P,F}。把式(46)代入式(44)～(45)得：

(47)

(48)

(49)

(50)

由于假定溫度保持不變，故式(47)～(50)中不含溫度變量。

在非飽和雙重孔隙介質研究中，為了簡化本構方程便于工程實用，通常假定裂隙和孔隙之間以及它們與各組分基質之間的力學性質相互獨立。此時由式(47)～(50)可得：

(51)

(52)

(53)

(54)

Pα=KRα?α

(55)

式中KRα=Kαα/φα0表示各組分基質的體積模量。假定非飽和雙重孔隙介質為各向同性材料，則KCC、KHH、KCrF和KCrP的具體表達式為：

(56)

(57)

(58)

(59)

式中I為二階單位張量，I4為四階單位張量，式中υC和EC分別是裂隙骨架的泊松比和彈性模量，υH和EH分別是孔隙骨架的泊松比和彈性模量，HCF和HrF分別為裂隙有效應力和吸力對裂隙飽和度的彈性模量，HHP和HrP分別為孔隙有效應力和吸力對孔隙飽和度的彈性模量。ζF和ζP按下式計算：

(60)

令KC=EC/[3(1-2υC)]和KH=EH/[3(1-2υH)]分別為裂隙和孔隙骨架的體積模量，Eb和υb為：

(61)

(62)

令Kb=Eb/[3(1-2υb)]。把式(56)～(59)代入到式(51)～(54)，對式(51)～(55)求逆后代入式(23)和式(35)～(38)，并把它們轉化為用σ和Pβγ表示的形式：

(63)

ζFL=-A12PT+(A22+φFL0/KRFL)PFL+A23PFG+

(64)

ζFG=-A13PT+A23PFL+(A33+φFG0/KRFG)PFG+

(65)

ζPL=-A14PT+A24PFL+A34PFG+(A44+

(66)

ζPG=-A15PT+A25PFL+A35PFG+A45PPL+

(67)

式中：

(68)

(69)

(70)

(71)

(72)

(73)

(74)

(75)

(76)

(77)

(78)

(79)

(80)

(81)

cFL、cFG、cPL和cPG根據文獻[8]的研究可表示為：

cFL=-cPL=ξL(PPL-PFL)

(82)

cFG=-cPG=ξG(PPG-PFG)

(83)

式(63)～(67)是小應變各向同性線彈性條件下非飽和雙重孔隙介質的固相本構方程、裂隙和孔隙中液氣兩相的滲入量本構方程。不難驗證：1)當非飽和雙重孔隙介質退化為飽和雙重孔隙介質時，φFG0、φPG0、sF、sP、1/HCF、1/HHP、1/HrP和1/HrF均為0，SrP0和SrF0等于1，把它們代入到式(63)～(64)和式(66)后可獲得Khalili線彈性方程；2)當忽略孔隙只考慮裂隙，非飽和雙重孔隙介質變為非飽和單重孔隙介質時，PPL、PPG、cFL、cPL、cFG、cPG、 1/EH、 1/KH和1/HrP均為0，把它們代入式(63)～(65)可得到非飽和單重孔隙介質線彈性方程[13]，這說明當雙重孔隙退化為單重孔隙或非飽和退化為飽和時，非飽和雙重孔隙介質的線彈性本構方程可退化為相應介質的線彈性本構方程。

把應力用應變表示時式(63)～(67)可變換得：

A13KbPFGI-A14KbPPLI-A15KbPPGI

(84)

(A23-A12A13Kb)PFG+(A24-A12A14Kb)PPL+

(85)

ζFG=A13KbεS∶I+(A23-A12A13Kb)PFL+

(A34-A13A14Kb)PPL+(A35-A13A15Kb)PPG-

(86)

ζPL=A14KbεS∶I+(A24-A12A14Kb)PFL+

(87)

ζPG=A15KbεS∶I+(A25-A12A15Kb)PFL+

(A35-A13A15Kb)PFG+(A45-A14A15Kb)PPL+

(88)

4 算例

4.1 基于本文理論的膨潤土試驗結果分析

在高放廢物深埋處置庫的概念設計中，需要采用膨潤土作緩沖材料充填于廢物罐和圍巖之間，阻止核素隨地下水遷移。經過研究比選，我國采用高廟子膨潤土作為核廢料處置庫的緩沖材料。由于膨潤土存在聚集體內的孔隙和顆粒間的裂隙，許多學者采用雙重孔隙介質來建立膨潤土的本構模型[6,8,18]。在這些模型中，大多根據孔隙水的性質，把孔隙視為被水飽和而裂隙被水部分飽和[6,8,18]，本文也遵循這一性質。

Morena等[18]假定當吸力大于20 MPa時膨潤土只有孔隙含水而裂隙不含水，總結了中國學者[19-20]研究高廟子膨潤土的系列試驗成果。筆者在Morena等[18]研究基礎之上，結合本文理論模型，重新整理了高廟子膨潤土的試驗數據。試驗時膨潤土固相、裂隙和孔隙體積分數為φS0=0.642、φF0=0.179和φP0=0.179，飽和孔隙骨架應變及其有效應力的關系見圖3，飽和裂隙骨架應變及其有效應力的關系見圖4，裂隙飽和度與膨脹應變之間的關系見圖5。

圖3 孔隙骨架豎向應變隨孔隙豎向有效應力變化Fig.3 Pore skeleton vertical strain changes with pore vertical effective stress

圖4 裂隙骨架豎向應變隨裂隙豎向有效應力變化Fig.4 Fracture skeleton vertical strain changes with fracture vertical effective stress

圖5 裂隙膨脹應變隨裂隙吸力變化Fig.5 Fracture swelling strain changes with fracture suction

裂隙土水特征曲線符合Van Genuchten模型：SrF=[1+(αSF)n]-m，式中m=0.825、n=1/(1-m)和α=2.16×10-4kPa-1。從圖3～5可以看出，膨潤土應力與應變存在較明顯的非線性關系。根據圖3、4和本文理論可得飽和排水條件下膨潤土應變隨總應力的變化曲線見圖6。圖6也給出根據單重孔隙介質理論獲得的膨潤土應變隨總應力曲線，從圖6可以看出本文提出的雙重孔隙介質理論比較符合試驗數據。

圖6 膨潤土豎向應變隨豎向總應力變化Fig.6 Bentonite vertical strain changes with vertical total stress

線彈性本構模型由于物理概念清楚，數學關系簡單，因此常用割線模量來近似建立土體的線彈性模型。泊松比按經驗值取為0.3，根據圖3～5提供的試驗數據和割線模型的計算方法，利用完全側限試驗的性質，膨潤土裂隙和孔隙骨架的力學參數見表1。

表1 非飽和雙重孔隙膨潤土MX80試樣的模型參數Tab.1 Model parameters of sample MX80 for unsaturated double-porosity bentonite

圖4給出的是飽和條件下獲得的試驗數據，非飽和條件下的楊氏模量比飽和條件下的有所提高，兩者可根據Alonso等[1]提出的計算公式進行換算[18]：EC(sF)=EC(0)/[0.836+0.164exp(-0.7sF)]，本算例裂隙吸力平均值約為1.2 MPa，故取吸力sF=1.2 MPa時的模量作為裂隙骨架的割線模量。試驗表明[21]高廟子膨潤土中水的滲透系數為5×10-11m/s，氣體的滲透系數為5×10-10m/s，由于孔隙假定飽和，無需KRPG值，故表中未給出相應值。根據表1和式(61)～(62)可得Eb=61 MPa，υb= 0.3，Kb=50.9 MPa。

4.2 核廢料緩沖材料的軸對稱固結控制方程

現在應用本文的線彈性本構方程來建立膨潤土作為核廢料儲藏罐的固結方程。把膨潤土緩沖設置簡化為圓柱軸對稱模型來分析，由于核廢料儲藏罐的體積較小，為簡便起見，把膨潤土視為實心圓柱體。注意到固結方程不考慮加速度和體力，首先把式(15)～ (16)按所有相相加得

divσ=0

(89)

在軸對稱條件下，由式(89)可得

(?σr/?r)+(σr-σθ)/r=0

(90)

(Kβγ/γw)Pβγ+nβγ0Wβγ=0

(91)

對式(91)求散度后應用ζβγ的定義并寫成軸對稱形式得

(92)

假定固相和水基質不可壓縮，孔隙被水飽和而裂隙被水部分飽和，把這些條件代入到式(84)～(88)后利用εr=?u/?r、εθ=u/r、式(90)和式(92)得：

(93)

(94)

(95)

(96)

式中θ=-(?u/?r+u/r)。設膨潤土半徑為R，由于線性系統的解答與土的初始應力無關，初始條件可等效地設置為：

ζβγ=0，Pβγ(R,t)=0，us(0,t)=0

(97)

(98)

式中p為作用在半徑R處圓柱體外側的外荷載，式(93)～(98)即是軸對稱條件下的固結控制方程和初邊值條件，如在推導式(93)～(98)方程中加上加速度項，就可以獲得波動控制方程。式(93)～(98)固結控制方程極其復雜，求解十分困難，但可以采用數值方法如差分法進行求解。

5 結 論

1)在考慮固相和流相基質變形的條件下，用嵌套思路推導了非飽和雙重孔隙介質的能量平衡方程。根據內能表達式中功共軛對的力學內涵揭示非飽和雙重孔隙介質本構模型的應變狀態變量為裂隙和孔隙骨架應變，裂隙和孔隙飽和度和各組分基質體應變；相應的應力狀態變量為裂隙和孔隙有效應力，裂隙和孔隙吸力和各組分基質壓力。

2)在小應變條件下，把固相應變分解為裂隙骨架應變、孔隙骨架應變和固相基質體應變之和。根據熱力學局部平衡假定，獲得非飽和孔隙介質的一般自由能勢函數本構方程。取自由能勢函數為狀態變量的二次多項式，獲得非飽和雙重孔隙介質的線彈性本構模型。然后根據混合物均勻化響應原理獲得各向同性線彈性本構方程。根據土工試驗獲得了非飽和膨潤土線彈性本構方程的模型參數。

3)把本文獲得的非飽和雙重孔隙介質的各向同性線彈性本構模型退化為飽和雙重孔隙介質和非飽和單重孔隙介質情形，可以獲得與前人相同的各向同性線彈性本構模型。把非飽和雙重孔隙介質線彈性本構方程與平衡方程和達西定理相結合，建立了非飽和膨潤土的軸對稱固結控制方程，可用于非飽和膨潤土防滲緩沖層的固結特性分析。