基于嵌套思路的饱和孔隙-裂隙介质本构理论

胡亚元

（浙江大学 滨海和城市岩土工程研究中心，浙江 杭州 310058）

自然界中，许多岩土材料具有两种不同尺度的孔隙，如裂隙黏土和岩体等.一种孔隙尺度比较小，通常仍称为孔隙，另一种孔隙尺度比较大，通常呈裂缝或扁平状，被称为裂隙.当孔隙和裂隙同时被一种流体占有时，就形成饱和孔隙-裂隙介质.近年来，随着水利水电、海底隧道、核废料储存以及海洋能源开发等工程大量建设，为了分析渗流和变形的流固耦合特性，饱和孔隙-裂隙介质的本构模型研究愈来愈受到工程力学界重视.Barenblatt 等[1]首先研究饱和孔隙-裂隙双重孔隙介质的本构特性.Khalili 等[2]、刘耀儒等[3]建立了各向同性饱和孔隙-裂隙介质的线弹性模型.蔡国庆等[4]和Zhao 等[5]建立了各向异性饱和孔隙-裂隙黏土的本构理论.张玉军等[6]创建了考虑裂隙产状等几何特性的孔隙-裂隙岩体的弹塑性模型.这些开创性成果有力地促进了饱和孔隙-裂隙介质力学本构理论的发展和应用.

在当前饱和孔隙-裂隙介质本构建模的研究文献中，针对同一个工程问题往往会创建出多种差异悬殊的本构模型.如何在各种模型中选择适合的饱和孔隙-裂隙介质本构模型成为工程师和学者首先遇到的难题.混合物理论从普适性的力学守恒定理出发研究孔隙-裂隙本构理论的普遍规律，具有严密的逻辑结构和明确的物理内涵，许多学者建议把混合物理论作为判定其他本构模型合理性的理论依据之一[7-11].Borja 等[7]和Zhang 等[8]根据混合物理论推导了饱和及非饱和孔隙-裂隙介质的能量平衡方程，并建立了饱和孔隙-裂隙介质线弹性本构模型，但该模型无法考虑裂隙与孔隙流相压力之差所导致的固相体积变化.Li 等[9-10]基于混合物理论推导了非饱和双孔隙膨胀土的外力功表达式，建立了非饱和双孔隙膨胀土的弹塑性本构模型；Guo 等[11]采用混合物理论建立了饱和及非饱和孔隙-裂隙介质的双有效应力弹塑性模型.然而，这些模型没有考虑固相和流相的材料变形，只适用于土体松散介质，无法适用于岩石和混凝土等非松散孔隙-裂隙介质[12-16].为了弥补上述缺陷，深刻揭示孔隙骨架应变和裂隙骨架应变在多孔介质流固耦合机制中的关键作用，便于利用均匀化响应原理相来建立相对简单实用的本构模型[14]，有必要对饱和孔隙-裂隙介质混合物理论作进一步深入研究.

鉴于此，笔者发现孔隙-裂隙介质可视为两个单重孔隙介质的嵌套叠加，即孔隙-裂隙介质可视为在单重裂隙介质的固相基质中嵌套了一个单重孔隙介质.本文从这一嵌套思路出发来研究饱和孔隙-裂隙介质的能量守恒方程和一般本构模型理论框架，从一般本构模型理论出发可推导饱和双重孔隙介质的线弹性方程，指导和校正当前饱和孔隙-裂隙介质的本构建模工作.

1 体积分数和密度

1.1 饱和孔隙-裂隙介质各组分体积分数和密度

饱和孔隙-裂隙介质是由固相、裂隙流相与孔隙流相组成的混合物.固相由S 表示，裂隙流相由F 表示，孔隙流相由P 表示.令α∈{S，F，P}为组分指征变量.φα为第α 组分的体积分数，ρα为第α 组分的平均密度，ρα为第α 组分的真实密度（或称材料密度），满足ρα=φαρα，则饱和孔隙-裂隙介质的总密度为ρ=ρS+ρF+ρP.根据体积分数的定义有：

1.2 基于嵌套思路的各组分体积分数和密度

本文把固相材料与孔隙流相组成的饱和单重孔隙介质称为饱和孔隙介质.当把饱和孔隙-裂隙介质中的固相材料和孔隙流相所构成的饱和孔隙介质视为一个整体时，此时只有裂隙被视为孔隙，本文把这种视角下的广义饱和单重孔隙介质称为饱和裂隙介质.这样，饱和孔隙-裂隙介质可看作在饱和裂隙介质的基质中嵌入饱和孔隙介质而成，而饱和孔隙-裂隙介质可视为两个单重孔隙介质的嵌套叠加.

根据上述嵌套思路，首先考虑饱和裂隙介质.饱和孔隙介质作为饱和裂隙介质的一个组分用SP 表示，它的体积分数为固相和孔隙流相体积分数之和φSP=φS+φP.根据式（1），在饱和裂隙介质中有：

然后考虑饱和孔隙介质.令β∈{S，P} 为饱和孔隙介质的组分指征变量，当饱和孔隙介质视为一个独立混合物时，则第β 组分在饱和孔隙介质中的体积分数为为饱和孔隙介质中固相的平均密度，则在饱和孔隙介质中有：

2 质量与动量守恒

2.1 质量守恒

令第α 组分的初始位置为Xα，t 时刻的空间位置为x，则每一组分的运动方程为x=xα（Xα，t），每一组分的速度和加速度可表示为：

对于定义在x 和t 上的标量场或矢量场Γα，基于α 组分的物质导数的定义为：

由于固相与裂隙中的流相、固相与孔隙中的流相均不存在质量交换，而裂隙中的流相和孔隙中的流相之间存在质量交换，则固相、裂隙流相与孔隙流相的质量守恒方程为：

式中：cF和cP分别表示裂隙流相与孔隙流相之间的质量交换率，满足cF+cP=0.

把固相作为饱和孔隙-裂隙介质混合物的参考构形，令裂隙流相和孔隙流相相对固相的扩散速度分别为WF=vF-vS和WP=vP-vS.把WF、WP、式（5）和ρα=φαρα代入式（6）～（8）得：

令div vr为饱和孔隙介质整体的体积变形率，与式（6）相类似，根据饱和孔隙介质中的质量守恒定律，有

将式（13）代入式（12）得：

2.2 动量和动量矩守恒

2.2.1 饱和孔隙-裂隙介质的Cauchy 应力张量

在饱和孔隙-裂隙介质混合物中，令σ 为混合物总Cauchy 应力张量，σα（α∈{S，F，P}）为第α 组分的Cauchy 应力张量，根据混合物理论有[7，16]

令总压力PT=-σ ∶I/3.固相材料真实压力PS与σS之间满足φSPS=-σS∶I/3，裂隙压力PF和孔隙孔压PP与其应力的关系满足φFPFI=-σF和φPPPI=-σP，利用上述关系和式（15）得：

根据嵌套思路，先考虑饱和裂隙介质.饱和孔隙介质作为饱和裂隙介质的组分，它的Cauchy 应力张量等于σSP=σS+σP.饱和裂隙介质总应力和各组分应力之间的关系由式（15）得：

令PSP=-σSP∶I（/3φSP）为饱和孔隙介质所受的真实压力，饱和裂隙介质总压力和各组分压力之间的关系由式（17）得：

再考虑饱和孔隙介质.如图1（c）和（d）所示，把饱和孔隙介质取为单元体，则固相和孔隙流相组分的Cauchy 应力张量分别为σrS=σS/φSP和σrP=σP/φSP，饱和孔隙介质混合物总应力张量为σr=σrS+σrP=σSP/φSP.故饱和孔隙介质在单元体上的总应力等于它在饱和裂隙介质中的真实应力，因而PSP=-σr∶I/3.饱和孔隙介质总压力和各组分压力之间的关系为：

图1 给出了饱和孔隙-裂隙介质总压力与各组分压力关系式.

图1 饱和孔隙-裂隙介质特征单元体示意图Fig.1 Schematic diagram for the representative volume element of saturated pore-fracture media

从图1 可以看出，图1（a）表示总压力PT作用在饱和孔隙-裂隙介质单元体上，图1（b）表示了式（16）和式（18）反映的饱和孔隙-裂隙介质和饱和裂隙介质的压力关系式，图1（c）表示PSP作用在饱和孔隙介质单元体上，图1（d）表示了式（19）反映的饱和孔隙介质的压力关系式.

2.2.2 饱和孔隙-裂隙介质的动量和动量矩守恒

工程界为便于应用，通常不考虑饱和孔隙-裂隙介质的微极介质特性，故第α 组分的动量矩供应量为0，利用固相、裂隙流相和孔隙流相的动量矩守恒方程可得应力张量σα（α∈{S，F，P}）是对称张量.

3 能量平衡方程

3.1 能量平衡方程

令qα、rα和分别为第α 组分的热流向量、外热供给量和能量供给量，ξα为第α 组分的内能密度，则固相、裂隙流相和孔隙流相的能量平衡方程为：

式中：DS=[grad vS+（grad vS）T]/2 为固相变形率；DF=[grad vF+（grad vF）T]/2 为裂隙流相变形率；DP=[grad vP+（grad vP）T]/2 为孔隙流相变形率.把式（23）～（25）相加可得：

从式（26）可知，等号右侧第一项的物理含义为各组分应变能变化率之和，利用σF=-φFPFI、σP=-φPPPI、WF和WP，各组分应变能变化率之和可得：

根据嵌套思路，先考虑饱和裂隙介质.令饱和裂隙介质Terzaghi 有效应力为=σ+PFI，利用式（2）、式（10）、式（13）和div vS=I ∶DS，把式（27）等式右边的前两项替换后可得：

式（32）表明DH的球应变速率与孔隙介质在裂隙介质中的体积分数相关.孔隙介质是裂隙介质的基质，在裂隙介质中起到骨架作用，因此本文把DH称为裂隙骨架变形率.式（34）表明DH的球应变速率与裂隙比改变率亦直接相关，即DH也可以采用裂隙介质的裂隙比来定义.同理，根据式（33）和式（35），DD可称为孔隙骨架变形率.把DD=Dr+（dSϑS/dt）I/3 和DH=DS-Dr代入式（31）后再把它代入到式（26）得：

在式（36）中，固相变形率DS被分为三部分：裂隙骨架变形率DH，孔隙骨架变形率DD和固相材料体应变率dSϑS/dt.在式（36）中，DH、DD、dSϑS/dt、dFϑF/dt 和dPϑP/dt 分别与和PP/ρF形成功共轭对.由热力学理论可知，在一般情况下，应选取裂隙骨架应变、孔隙骨架应变、固相材料体应变、裂隙流相材料体应变和孔隙流相材料体应变作为饱和孔隙-裂隙介质本构模型的应变状态变量；选取单位密度上的裂隙介质有效应力、孔隙介质有效应力、固相材料真实压力、裂隙孔压和孔隙孔压作为应力状态变量.

3.2 混合物均匀化响应原理

为了适应工程应用，工程界常常利用混合物均匀化响应原理来简化混合物的本构关系.混合物均匀化响应原理的内容为[14]：当混合物单元体承受外荷载时，若混合物单元体中每一点的真实应变增量（或速率）相等，则该混合物单元体等效于单相均匀单元体，即单元体内每一点处的真实应力增量（或加荷速率）也相等；反之也然.在Khalili 等[2]、陈正汉[17]、陈勉和陈至达[18]推导各种饱和和非饱和混合物本构关系时，混合物均匀化响应原理曾发挥了至关重要的作用.

现在应用混合物均匀化原理对能量平衡方程的合理形式作分析.根据嵌套思路，先应用混合物均匀化原理分析饱和裂隙介质.令裂隙介质有效压力为∶I/3.在保持总压力PT与裂隙孔压PF增速相等的情况下，饱和裂隙介质每一点的真实应力增速相等，此时裂隙介质有效压力增速为零，但在饱和孔隙介质变形不可忽略的情况下，饱和孔隙介质将产生大小为Dr∶I≠0 的体应变速率.根据混合物均匀化响应原理，饱和裂隙介质每一点应力加荷速率相等时，其应变速率也相等，则饱和裂隙介质的体应变速率为DS∶I=Dr∶I ≠0，故当时DS≠0，这意味着DS不仅与有关还与饱和裂隙介质的其他应力有关，故选取DS作为的功共轭变量不利于饱和裂隙介质的本构建模工作.因为，当dt=0 时，DH=DS－Dr=0，所以本文选取DH作为的功共轭变量，这样当混合物均匀化原理成立时，DH只与有关而与饱和裂隙介质的其他应力无关，从而根据弹性互易定理可以推断出DH与其他应变相互解耦的结论.

再应用混合物均匀化响应原理分析饱和孔隙介质.与饱和裂隙介质的分析相类似，当混合物均匀化响应原理成立时，Dr不仅与有关还与饱和孔隙介质的其他应力有关，故选取Dr作为的功共轭变量同样不利于本构建模工作.由于当时，DD=Dr+（dSϑS/dt）I/3=0，故选取DD作为的共轭变量.当均匀化响应原理成立时，同理可以得出Dr与其他应变相互解耦的结论，这也是本文要把应变能写为式（31）和能量平衡方程写为式（36）的原因.

3.3 能量平衡方程的退化

式（38）与饱和多孔介质的能量守恒方程完全一致[14].

4 一般势函数本构方程

4.1 有限应变情况下的一般本构方程

在有限应变情况下，利用连续介质力学[19]中变形梯度的分解方法，先将裂隙介质变形梯度分解为孔隙介质变形梯度与裂隙骨架变形梯度的乘积；后将孔隙介质变形梯度分解为固相材料变形梯度与孔隙骨架变形梯度的乘积，如图2 所示，图中F 表示变形梯度.

图2 孔隙-裂隙介质变形梯度示意图Fig.2 Schematic diagram for the deformation gradient of pore-fracture media

为了工程应用，通常做如下简化：孔隙-裂隙介质的固相内能、裂隙流相材料内能和孔隙流相材料内能之间相互独立，即ξS=ξS（UH，ED，ϑS，ηS），ξF=ξF（ϑF，ηF）和ξP=ξP（ϑP，ηP），ηα为各组分熵密度，并假定固相材料、裂隙流相材料和孔隙流相材料具有相同温度θ.则根据热力学局部平衡条件以及UH、ED、ϑS、ϑF、和ϑP相互独立的性质，有

式（40）和式（41）便是有限应变情况下饱和孔隙-裂隙介质的内能势函数一般本构方程.由于内能是一种自由能，反映的是弹性性质，故式（40）和式（41）为饱和孔隙-裂隙介质的一般弹性本构方程.把式（40）和（41）代入式（39）得：

将式（42）与非平衡态热力学相结合，按照文献[14]推导过程可进一步获得饱和孔隙-裂隙介质的塑性本构方程，受篇幅限制本文不再赘述.

4.2 小应变情况下的一般本构方程

4.2.1 小应变情况下各组分应变计算公式

在小应变情况下可略去高次项，令εS为固相应变张量ES的近似值，εH为裂隙骨架应变张量UH的近似值[19]，εr为孔隙介质应变张量Er的近似值，εD为孔隙骨架应变张量ED的近似值，此时得[19]：

式中：uH=xS-xr，ur=xr-XS，uD=xr-xM，xM是材料变形后的位置.注意到ϑS以压为正，ε 以拉为正.对比式（43）～式（47）可得：

由式（48）可知，在小应变情况下，饱和孔隙-裂隙介质的固相应变εS可以分解为裂隙骨架应变εH、孔隙骨架应变εD与固相材料体应变ϑS之和.

在小应变情况下，令εSV、εHV和εDV分别为固相、裂隙骨架和孔隙骨架体应变，对式（48）取迹得：

根据式（34）和式（35），可得εHV和εDV的表达式分别为：

式中：φSP0为饱和孔隙介质整体的初始体积分数：为饱和孔隙介质中固相组分初始体积分数.

接下来推导流相体应变计算式.首先，在小应变情况下，ϑF和ϑP可近似地简化为：

式（58）和（59）便是小应变情况下裂隙流相体应变和孔隙流相体应变的计算式.从式（49）、式（58）和（59）可以明显地看出，裂隙骨架体应变同时影响固相体应变、裂隙流相体应变和孔隙流相体应变，孔隙骨架体应变同时影响固相体应变和孔隙流相体应变.故固相、孔隙流相和裂隙流相之间必然存在受力变形耦合，它们是通过裂隙骨架和孔隙骨架应变进行传递和协同的.

4.2.2 一般本构方程

令ξ*=ρS0ξS+ρF0ξF+ρP0ξP，式中：ρS0为固相初始平均密度；ρF0为裂隙流相初始平均密度；ρP0为孔隙流相初始平均密度.小应变条件下可略去高次项，有

根据热力学局部平衡假定，式（61）中的内能可表示为ξ*=ξ*（η，εH，εD，ϑα），对它求全微分得：

根据热力学局部平衡假定，对比式（61）和式（62）得：

式（63）便是饱和孔隙-裂隙介质在小应变情况下的一般弹性本构方程.将式（64）与非平衡态热力学相结合，按照文献[14]的推导过程可获得饱和孔隙-裂隙介质的塑性本构方程，受篇幅限制，本文不再赘述.

4.3 一般本构方程的退化和验证

4.3.1 小应变各向同性线弹性本构方程

引入Helmhotlz 自由能ψ*（θ，εH，εD，ϑα），它等于ψ*（θ，εH，εD，ϑα）=ξ*（η，εH，εD，ϑα）-θη，对它求全微分后把式（63）代入得：

令饱和孔隙-裂隙介质的初始平衡状态为（θ，εH，εD，ϑα）=（θ0，0，0，0），在受到微小扰动后，到达一个新的平衡状态（θ0+θΔ，εH，εD，ϑα），θΔ为温度增量.根据势函数本构方程的一般性质可知，若自由能函数取为状态变量的二次多项式函数，可获得线弹性本构关系.故Helmhotlz 自由能取为：

式 中：KHH、KDD、Kαχ、Kθθ、KHD、KHα、KDα、KHθ、KDθ、Kθα为模型的弹性系数；Kαχ=Kχα，α∈{S，F，P}，χ∈{S，F，P}.把式（66）代入式（65）得：

式中：KDH为KHD的转置.

在当前饱和孔隙-裂隙介质研究中，绝大多数研究均假定：1）温度不变，即θΔ=0；2）裂隙与孔隙中流相材料的本构关系与其单独存在时的本构关系相同；3）孔隙骨架和裂隙骨架存在变形耦合；4）固相介质为各向同性材料.由假定2）可知固相变形与流相材料变形相互解耦.根据假定4），KHH和KDD为各向同性张量.根据上述4 个假定，式（67）～式（69）可分别表示为：

式中：I 为二阶单位张量；I4为四阶单位张量；νHH为裂隙骨架自身的泊松比；EHH为裂隙骨架自身的弹性模量；νHD为孔隙和裂隙骨架的耦合泊松比；EHD为孔隙和裂隙骨架的耦合弹性模量；νDD为裂隙骨架的泊松比；EDD为孔隙骨架的弹性模量；KRα、α∈{S，F，P}分别为固相材料和流相材料的体积模量.AV和AS分别为：

对式（71）和（72）求逆后得：

式（84）即为固相线弹性本构方程.对式（84）求迹并利用PT=-σ ∶I/3得固相体应变为：

在实际工程中，比较关心的是流体流出或流入孔隙-裂隙介质的流量，定义裂隙流相渗入量为ζF=φF0（εSV-εFV），孔隙流相渗入量为ζP=φP0（εSV-εPV）[20]，利用式（49）、式（58）和（59），ζF和ζP的表达式分别为：

式（88）和（89）分别为本文获得的裂隙流相和孔隙流相渗入量本构方程.

4.3.2 与Khalili 方程的对比与验证

Khalili 等[2]假定1/KHD=0，获得的固相本构方程、裂隙流相渗流方程和孔隙流相渗流方程分别为：

式中：kF和kP分别为裂隙骨架与孔隙骨架的各向同性渗透系数；μ为流相材料的黏滞系数：为拉普拉斯算子.对式（90）求εS关于σ 的反函数可得：

注意到Khalili 等[2]推导式（91）和（92）时假定裂隙和孔隙流相渗流满足达西定理.根据达西定理可得：

把式（94）和（95）代入式（91）和（92）后求时间积分得：

显然，当1/KHD=0 时，式（84）、式（88）和（89）与式（93）、式（96）和（97）完全一致，说明从本文的自由能势函数一般本构方程出发可以获得与Khalili 等相同的线弹性本构模型.Khalili 等把他们的线弹性本构模型用于裂隙黏土的固结分析，获得了与试验数据相一致的理论分析结果[2，15].这说明从本文的一般本构方程出发可获得经过试验验证的本构模型.

5 结论

1）在考虑固相和流相材料变形的条件下，以嵌套思路推导了饱和孔隙-裂隙介质的能量平衡方程.确定了饱和孔隙-裂隙介质本构方程的应变状态变量是裂隙骨架应变、孔隙骨架应变、固相材料体应变、裂隙流相材料体应变和孔隙流相材料体应变；应力状态变量是单位密度上的裂隙介质有效应力、孔隙介质有效应力、固相材料真实压力、裂隙孔压和孔隙孔压.

2）在小应变情况下，固相应变可分解为裂隙骨架应变、孔隙骨架应变和固相材料体应变之和.获得有限应变和小应变条件下的饱和孔隙-裂隙介质的自由能势函数一般本构方程.

3）当混合物均匀化响应原理成立时，裂隙骨架、孔隙骨架和固相材料的本构模型相互解耦；当裂隙与孔隙中流相材料的本构关系与纯流相本构关系相同时，固相与流相材料变形相互解耦.当上述两个性质均成立时，裂隙骨架应变唯一决定裂隙介质有效应力、孔隙骨架应变唯一决定孔隙介质有效应力、固相材料体应变唯一决定固相材料真实压力、裂隙流相材料体应变唯一决定裂隙孔压和孔隙流相材料体应变唯一决定孔隙孔压.运用这些本构性质可以简化本构关系的复杂程度，有利于工程应用.

4）当自由能势函数取为状态变量的二次多项式时，获得孔隙骨架和裂隙骨架相互耦合的各相同性线弹性本构方程，当孔隙骨架和裂隙骨架变形解耦时，该线弹性方程退化为饱和孔隙-裂隙介质Khalili线弹性方程.Khalili 等利用他们提出的线弹性本构方程获得与试验数据相一致的理论分析结果[2，15]，这说明本文基于一般势函数的本构方程理论框架可以指导饱和孔隙-裂隙介质的具体本构建模工作.