双应力变量的饱和介质非线性体积本构关系

根据工程影响范围[1]的大小，岩体裂隙可分为众多网格状、宽度小、间距近的低序次微小裂隙，以及数目不多的宽度大、间距远的高序次大中型裂缝和断层.大部分学者把诸如砂岩、泥岩和黏土岩等包含大量微小裂隙且孔隙率较大的饱和岩体视为多孔介质，并采用拟连续介质力学研究其本构关系和力学特性[2-9].由文献[10-13]可知，饱和岩体力学需要考虑固相基质变形(在饱和岩体力学中称为岩体矿物变形)的性质，而目前这一性质很难体现在饱和岩体非线性本构模型中，因此限制了饱和岩体多场耦合理论的应用.为突破这一限制，亟需建立更加合理和完善的饱和多孔介质本构理论体系.

选择合适的状态变量可以有效降低饱和多孔介质本构方程的建模难度.例如，Terzaghi[14]利用饱和土体固相基质变形(在土力学中称为土颗粒变形)可忽略不计的性质提出Terzaghi有效应力原理，并以此作为饱和土体本构方程的唯一状态变量，成功创建了饱和土力学.鉴于Terzaghi有效应力原理在饱和土力学中的成功应用，岩土工程界试图沿用单一的有效应力变量来建立饱和多孔介质本构理论.Biot[15]指出当饱和多孔介质(如饱和岩体或混凝土)的固相基质变形不能忽略时，需要对Terzaghi有效应力公式中的孔压进行折减才能适用于一般饱和多孔介质，该折减系数称为Biot系数.Skempton[16]分别将摩擦角、剪切抗力角、面积比、固相基质体积模量和固体骨架体积模量等力学参数引入Biot系数公式，并针对固相变形和强度等效原则提出了两种不同的有效应力表达式.然而，现代本构理论认为强度问题是变形问题的延伸，它们采用的应力状态变量应当具有统一性，即应力状态变量不会随着强度与变形问题的变化而发生改变.因此，Skempton提出“针对变形和强度问题需要采用不同的有效应力表达式”的研究结论暗示饱和多孔介质的全部力学性质不能只通过一个有效应力变量来描述.Borja[12]为了获得单应力变量的饱和多孔介质力学理论，提出非线性的Biot系数公式.但是，Borja 在推导Biot系数公式过程中，假定固相(骨架)体应变与有效应力无关，这个假设在饱和土中与Terzaghi有效应力原理相违背.其他岩土学者也针对具体岩土问题提出了形式不一的有效应力表达式[13]，虽然这些表达式适应了特定的使用范围，但其普适性局限问题接踵而来.饱和岩体试验表明，采用单一有效应力变量很难表征岩体复杂多变的非线性性质.混合物理论公设饱和多孔介质的力学性质由固流两相共同决定，故主张采用双应力变量来研究饱和多孔介质的本构模型.如：De Boer[17]采用总应力和孔压建立了考虑体积分数变化的饱和弹塑性多孔介质的混合物理论；Schrefler[18]建议把体积平均化理论引入混合物理论来定义固流两相的应力和应变，研究包括饱和岩体在内的饱和多相孔隙介质；Borja[12]基于平均化理论，把宏细观力学与混合物理论相结合，采用有效应力和孔压获得了饱和多孔介质的能量守恒方程和机械功公式.这些研究结果表明：对于必须考虑基质变形的饱和多孔介质，采用双应力或双应变状态变量可更加灵活地模拟饱和多孔介质的力学性质，解释固流两相的相互作用机理，从而建立更加完善的本构关系理论体系.

本文基于工程混合理论[19]推导和分析具有普适性质的饱和多孔介质体积变形功守恒方程；提出固相基质压力与Terzaghi有效应力共同决定饱和多孔介质固相体应变的双应力原理；结合模型试验，建立基于双应力变量的饱和多孔介质非线性体积本构方程.在此基础上，分析Biot系数、孔压系数、孔隙比等随Terzaghi有效应力的变化规律，开辟了研究基质变形的非线性饱和多孔介质本构性质的新方法.

1 体应变本构方程的建模理论

根据经典饱和多孔介质混合物理论[12,20]，将固流两相按照体积分数平均化后得到的连续体构形产生的应变称为组分应变，将固流两相实际存在的细观真实构形产生的应变称为组分基质应变.假定饱和多孔介质的固流两相之间不存在质量交换.用S表示固相，F表示流相，对于饱和多孔介质，体积分数满足：

φS+φF=1

(1)

图1 饱和多孔介质表征单元体图
Fig.1 The diagram for the representative volume element of saturated porous media

设σSm和σFm分别为混合物理论[19]中固流两相承受的Cauchy球应力；pS为固相基质承受的球应力，简称固相基质压力，pS=σSm/φS；pF为流相基质承受的球应力，简称流相基质压力或孔压，pF=σFm/φF；σm为饱和多孔介质承受的总球应力.单元体受力情况如图1所示.针对选取的横截面ABCDE，根据混合物理论有[19]：

σm=σSm+σFm=φSpS+φFpF

(2)

设初始时刻(t0)固相和流相所占的真实体积分别为VRS0和VRF0，体积分数分别为φS0和φF0，则混合物的总体积为V0=VRS0+VRF0.按体积分数平均化后固相体积为VS0，VS0=VRS0/φS0，流相体积为VF0，VF0=VRF0/φF0.t时刻，固流两相在各自应力的作用下体积分别变为VS和VF，固流两相基质体积分别变为VRS和VRF.由于流相在混合物中存在渗流运动，所以t时刻VS一般不等于VF.

设：初始固相基质真实密度为ρRS0，初始固相组分的相密度为ρS0，ρS0=φS0ρRS0；初始流相基质真实密度为ρRF0，初始流相组分的相密度为ρF0，ρF0=φF0ρRF0；参数ρα(α=RS0,S0,RF0,F0)在Cauchy球应力作用下变化为ρβ(β=RS,S,RF,F).混合物理论是按平均化后的组分体积定义组分体应变的，令固相和流相的体应变分别为εSV和εFV，由于变形前后组分的质量保持不变，有[19]

εSV=-ln(VS/VS0)=ln(ρS/ρS0)

(3)

εFV=-ln(VF/VF0)=ln(ρF/ρF0)

(4)

当以压为正时，根据连续介质力学理论，固流两相基质体应变εRSV、εRFV的定义为[21]

εRSV=-ln(VRS/VRS0)=ln(ρRS/ρRS0)

(5)

εRFV=-ln(VRF/VRF0)=ln(ρRF/ρRF0)

(6)

把由固相体积分数变化引起的固相体应变定义为固相体积分数应变εSf[22-23]，

(7)

式中：e和e0分别为多孔固相的当前和初始孔隙比，且有

φS=1/(1+e)， φF=e/(1+e)

由式(3)～(7)可得

εSV=εSf+εRSV

(8)

εFV=ln(e/e0)+εSf+εRFV

(9)

根据式(5)、(7)和(8)，固相体应变可分解为固相基质体应变和固相体积分数应变之和，如图2所示.

图2 固相体应变分解图
Fig.2 Decomposition diagram of the volume strain of solid phase

根据混合物理论[12,19]，饱和多孔介质体积变形功率为固、流相球应力与体应变率乘积之和：

(10)

将式(8)和(9)对时间微分后代入式(10)，并联立式(1)和(2)可得

(11)

式(11)表明：固流两相所受的基质压力与基质体应变为功共轭对；Terzaghi有效应力与固相体积分数应变为功共轭对；饱和多孔介质混合物的总体积变形功等于固流两相基质体积变形功与固相体积分数变形功之和.

假定饱和多孔介质单位质量的变形功由3部分组成：① 固相单位密度内所含的固相细观变形引起的变形功WRS，只与εRSV有关；② 流相单位质量内所含的流相细观变形引起的变形功WRF，只与εRFV有关；③ 固相单位密度内所含的固相体积分数变化引起的变形功WSf，只与εSf有关.根据混合物理论[19]，有

1.4.6 归一化积雪指数(NDSI) 归一化积雪指数(NDSI)是提取积雪信息的一种有效方法，其算法较合理，分类精度高，具有普遍的操作意义。NDSI类似于归一化植被指数(NDVI)，对大范围的光照条件不敏感，对大气作用可使其局地归一化并且不依赖于单通道的反射，NDSI的计算公式：

(12)

pS=ρRS0exp(εRSV)(dWRS/dεRSV)

(13)

pF=ρRF0exp(εRFV)(dWRF/dεRFV)

(14)

(15)

在小应变条件下，有φS≈φS0和φF≈φF0.根据式(2)和式(8)、(9)有

σm=φS0pS+φF0pF

(16)

εSV=εSf+εRSV

(17)

εFV=-φS0εSf/φF0+εRFV

(18)

式中:

εSf=(e0-e)/(1+e0)

(19)

εRSV=(ρRS-ρRS0)/ρRS0

(20)

εRFV=(ρRF-ρRF0)/ρRF0

(21)

式(19)与不考虑固相基质压缩性(即不考虑土颗粒压缩性)的饱和土固相变形公式[24]一致.为了简洁并参照一般教科书[21]在小应变条件下的约定，在式(16)～(21)中把“≈”写成“=”.在小应变条件下，εRSV、εRFV和εSf均可视为无穷小量，ρRS0、ρRF0和ρS0是常量，因此由式(13)～(15)可得

pS=d(ρRS0WRS)/dεRSV

(22)

pF=d(ρRF0WRF)/dεRFV

(23)

(24)

式(22)～(24)表明：在小应变条件下，固相基质体应变唯一决定固相基质压力；流相基质体应变唯一决定流相基质压力；固相体积分数应变唯一决定Terzaghi有效应力，反之亦然.式(17)表明固相体应变等于固相体积分数应变和固相基质体应变之和.因此，固相体应变由Terzaghi有效应力和固相基质压力共同决定，即双应力原理.若忽略固流两相基质的压缩变形特性，则有εRSV =εRFV =0，εSV=εSf，此时双应力原理退化为有效应力原理.

综上所述，一个完整的饱和多孔介质本构方程除需确定固相应变外还需确定流相应变.由式(17)、(18)和式(22)～(24)可知，需同时利用固相基质体应变、固相体积分数应变和流相基质体应变(或固相基质压力、Terzaghi有效应力和孔隙流体压力)，才能反映小应变条件下饱和多孔介质的全部受力变形性质.

2 试验模拟与分析

2.1 饱和多孔介质双应力-应变关系拟合

鉴于岩土工程界对于饱和多孔介质有效应力公式存在诸多分歧的问题，Lade等[13]开展了等向压缩模型试验，以期推导正确的有效应力表达式.他们将强压缩性材料胶结成大孔隙立方体固相骨架后完全浸水形成饱和多孔介质模型.在等向压缩条件下通过施加不同的总球应力和孔压，测量了饱和多孔介质模型的固相压缩体积变形和固相基质压缩体积变形.本文根据双应力原理，结合试验数据建立以Terzaghi有效应力与固相基质压力为状态变量的饱和多孔介质本构模型.

试样初始总体积为 203.2 cm3，φF0=0.436[13].在小应变条件下，(φF-φF0)相对于φF0为无穷小量，可以忽略不计，因此φF=φF0+(φF-φF0)≈φF0.将原始测量数据以及各基质应力与对应的基质应变数据绘制于图3.通过式(22)所示性质进行曲线拟合，结果表明εRSV与pS两者之间呈线性关系：

εRSV=CSpS

(25)

式中：CS为固相基质的压缩系数，CS=2.2×10-5 kPa-1.拟合度R2=0.922.

由式(24)可知，在小应变情况下εSf仅与相关.将εSf随变化的试验数据绘制于图4中并通过曲线拟合发现两者呈双曲线关系：

(26)

式中：a和b为模型参数.根据图4中的数据进行回归分析可得a=83，b=9 136，R2=0.96.

图3 固相基质压力与对应基质应变的关系曲线
Fig.3 Relationship curve between solid matrix strain and solid matrix pressure

图4 固相体积分数应变与Terzaghi有效应力关系曲线
Fig.4 Relationship curve between solid volume fraction strain and Terzaghi effective stress

根据式(17)、(25)和(26)可得εSV为

(27)

式中：式(27)是通过双应力变量获得的固相体应变本构关系式.根据式(25)～(27)，只要通过试验获得CS、a和b值，就能够通过Terzaghi有效应力和固相基质压力计算出固相体积分数应变、固相基质体应变和固相体应变.图5为本试验固相体应变与Terzaghi有效应力以及固相基质压力的关系曲面.由图5可见，固相体应变随着双应力变量的增加逐渐递增.

图5 固相体应变随双应力变化曲面
Fig.5 The variation surface of solid phase strain with double stresses

2.2 基于双应力变量模型的Biot系数

Biot[15]认为当考虑固相基质压缩性时，只要引入Biot系数(η)对Terzaghi有效应力公式[14]中的孔压进行修正，就仍可以采用一个有效应力表征一般饱和多孔介质的固相变形特性.Skempton[16]根据固相变形等价原则获得了线弹性多孔介质的Biot系数表达式.本文借鉴Skempton方法，根据双应力变量建立的饱和多孔介质本构模型推导非线性条件下的Biot系数表达式.

按Biot理论，固相体应变可以采用修正后的有效应力表示为

εSV=Cb(σm-η pF)

(28)

式中：Cb为固相骨架压缩系数.根据修正后有效应力表征的固相体应变与双应力表征的固相体应变等效的性质，联立式(27)和(28)得

(29)

由于式(29)对任意的σm和pF均成立，所以

把本模型测定的CS、a和b的值代入式(30)和(31)，可以得到Cb与η随Terzaghi有效应力的变化关系，如图6所示.

图6 Cb和η与的关系曲线
Fig.6 Relationship curve between Cb, η and 

由图6可见，Cb与η均随的增大而减小.出现此现象的原因主要有两方面.首先，随着外荷载增加，固相骨架逐渐压缩紧密，孔隙比减小，产生相同压缩应变所需荷载更大，导致Cb减小.其次，随着外荷载增加，孔隙流体逐渐从骨架排出，孔隙压力对等效有效应力影响所占的权重逐渐减小，故表征该权重的η随外荷载的增加而逐渐减小.因此，对于考虑固相基质压缩性的非线性饱和多孔介质，外荷载作用下孔隙逐渐被压缩时，其Cb与η不再是常数，而是呈现出随 增加而逐渐减小的变化趋势.

2.3 基于双应力变量模型的孔压系数

在岩土力学中把不排水三轴试验条件下孔压增量与总压增量之间的比值定义为孔压系数(B)，现采用双应力变量模型来研究B的表达式.

假设某一时刻单元体体积为V，小应变条件下φS≈φS0，φF≈φF0.在不排水的条件下施加一个外荷载，使总球应力增量为Δσm，产生的孔隙水压力增量为ΔpF，固相基质压力的增量为ΔpS，ΔpS=(Δσm-φF0ΔpF)/(1-φF0)，Terzaghi有效应力的增量为孔隙水(流相基质)压缩系数取CF，CF=4.8×10-7 kPa-1 [13].根据不排水条件下的体积变形关系可得固相体积压缩量等于固流两相基质体积压缩量之和[25]：

ΔεSVV=ΔεRSV(φS0V)+ΔεRFV(φF0V)

(32)

根据式(27)，双应力表征的固相体应变增量为

ΔεSV=

(33)

固流两相基质的体应变为

ΔεRSV=CSΔpS， ΔεRFV=CFΔpF

(34)

把式(33)和(34)代入到式(32)得

(35)

式中:从式(35)可以看出:对于非线性饱和多孔介质，B不但与有关，而且与有效应力增量有关.可见B的变化比较复杂.为了简化B的变化规律，本文提出切线孔压系数(Bt)的概念，并定义由Bt的定义和式(35)得

(36)

由式(36)可知，当CF<CS时，Bt>1；当CF>CS时，Bt<1.

在实际工程中，饱和多孔介质的孔隙液体中溶解了少量气体，一般把饱和度Sr∈(0.95,1)的饱和多孔介质称为准饱和多孔介质.如果反压后试验时间过长，模型试验的饱和试样并不是完全饱和的，而是处于准饱和状态(0.95<Sr<1) [26].因此，考虑准饱和状态下孔压系数变化规律具有实际工程意义[26-27].对于准饱和多孔介质，孔隙流体包含了液体与气体，其压缩系数简化为[28]

CF=Cg(1-Sr)+ClSr

(37)

式中：Cg为气体压缩系数，Cg=4.94×10-3 kPa-1 [25]；Cl为液体水压缩系数，Cl=4.8×10-7 kPa-1 [13,28].

对于本模型试验，Lade等[13]采用轻木、椴木、聚丙烯3种材料模拟饱和多孔介质的固相基质，在完全饱和状态下3种材料的CS均大于CF.将本模型试验测定的a和b值以及式(37)代入式(36)，得到准饱和状态下饱和度对孔压系数的影响规律，如图7所示.可以看出：Sr>99.7% 时Bt均大于 1.0.由于孔隙中含有微量气体，导致孔隙流体压缩系数偏大.当Sr=99.7% 时，孔隙流体压缩系数与固相基质压缩系数非常接近，此时Bt接近于 1.0.随着Sr进一步降低，含气流体的压缩系数大于固相基质压缩系数，此时，Bt小于1.0.这一结果与式(35)和(36)显示的变化规律一致，说明孔隙流体中含气量的微小变化会显著改变Bt的变化规律.

图7 切线孔压系数与Terzaghi有效应力关系
Fig.7 Relationship curve between tangent pore pressure coefficient and Terzaghi effective stress

2.4 基于双应力变量模型的孔隙比和孔隙率变化规律

根据双应力原理可以得到固相体积分数应变，从而推导出饱和多孔介质的孔隙比和孔隙率，对于研究饱和多孔介质变形过程中的孔隙变化具有重要意义.在小应变条件下根据式(19)和(26)以及e0=φF0/(1-φF0)，可以获得e和φF随Terzaghi有效应力的变化公式为

(38)

(39)

把本模型试验测定的a和b值代入式(38)和(39)，得到e与φF随着Terzaghi有效应力变化的关系，如图8所示.从图中可以看出：随着外荷载增加，固相骨架逐渐压缩紧密，孔隙逐渐减少；减小，导致骨架变形逐渐减少，孔隙比减小变缓并逐渐趋于稳定.由图中还可知，孔隙变化量较小(孔隙比从 0.773 减小到 0.754，孔隙率从 0.436 减小到 0.430)，属于小变形范围，此时有e≈e0和φS≈φS0.

图8 e和φF与关系曲线
Fig.8 Relationship curve between e, φF and 

3 结论

本文基于混合物理论以及连续介质力学，根据不同研究尺度范畴，将固流两相在饱和多孔介质中的存在构形分为细观真实构形与连续体构形，并证明了固相体应变等于固相基质体应变和固相体积分数应变之和.获得的研究结果如下.

(1) 提出小应变条件下考虑基质压缩性的饱和多孔介质力学的双应力原理：Terzaghi有效应力唯一决定固相体积分数应变；固相基质压力唯一决定固相基质体应变；Terzaghi有效应力与固相基质压力共同决定固相体应变.当忽略固相基质压缩性时，固相基质体应变为0，双应力原理退化为有效应力原理.

(2) 根据双应力原理，结合模型试验数据进行拟合，得到固相体积分数应变和Terzaghi有效应力之间以及固相基质体应变和固相基质压力之间的关系方程，确定了模型参数，从而建立了饱和多孔介质的双应力变量本构模型.

(3) 基于双应力变量本构模型，得出饱和多孔介质固相骨架压缩系数、Biot系数、孔压系数，以及孔隙比和孔隙率的表达式.其中固相骨架压缩系数与Biot系数随Terzaghi有效应力的增加呈双曲线递减变化规律；孔压系数与固相基质压缩系数以及多孔介质的准饱和度有关；随着外荷载增加，孔隙逐渐减小，固相骨架压缩系数亦随之减小，孔隙比与孔隙率减小速度减慢并呈现出逐渐稳定趋势.