非饱和双重孔隙介质的势函数本构方程及应用

胡亚元

(浙江大学 滨海和城市岩土工程研究中心，杭州 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)把本文获得的非饱和双重孔隙介质的各向同性线弹性本构模型退化为饱和双重孔隙介质和非饱和单重孔隙介质情形，可以获得与前人相同的各向同性线弹性本构模型。把非饱和双重孔隙介质线弹性本构方程与平衡方程和达西定理相结合，建立了非饱和膨润土的轴对称固结控制方程，可用于非饱和膨润土防渗缓冲层的固结特性分析。