坝基的渗流变形和稳定一直是关键的工程问题^[1-3]. 坝底基岩包含众多网格状的低序次微小孔隙和数目不多的高序次裂隙，具有典型的孔隙-裂隙结构特征. 孔隙不易压缩，对压力不敏感，渗透能力弱，但具有较强的蓄水能力，在岩体中主要起贮存作用；裂隙易压缩，压敏性较高，渗透能力较强，在岩体中充当渗流通道作用^[4-7]. 孔隙和裂隙在尺寸、变形性质和渗透特性中存在显著差异，可根据岩体结构面的分类标准，将孔隙系统归类为Ⅴ级结构面，将裂隙归类为Ⅳ级结构面^[8]. 针对上述工程岩体的多尺度特征，Barrenblatt等^[4]提出双重孔隙介质概念，在固相刚性条件下建立饱和孔隙-裂隙岩体的渗流耦合方程. Elsworth等^[6]应用孔隙和裂隙双有效应力原理直观地反映了孔、裂隙骨架的变形，但没有明确孔隙和裂隙骨架应变的数学定义，根据试验确定的模型参数不符合模型的力学含义. Berryman等^[7]唯象地建立了全耦合的线弹性孔隙-裂隙本构模型，该模型参数众多，在实际工程中运用困难. 张玉军^[9]在Elsworth等^[6]模型的基础上考虑了温度影响，建立孔隙-裂隙岩体多场耦合模型，但没有给出模型参数的取值依据. 上述模型在创建过程中往往基于直觉将单孔介质弹性理论推广到孔隙-裂隙介质中，在建立固相总应变与固相材料应变、孔隙骨架应变和裂隙骨架应变之间的关系时缺少严密的数学推导.

很多学者采用混合物理论来研究饱和孔隙-裂隙介质流固耦合的变形规律，Wilson等^[5,10-11]基于混合物理论，直接以组分应力或应变作为状态变量，建立考虑固相变形的流固耦合方程. 这些本构方程的模型参数在工程和试验中难以直接测量，限制了它们在工程中的应用. Borja等^[12]基于混合物理论建立孔隙-裂隙非饱和介质的一般本构理论. 在不考虑固流两相材料变形的条件下，Li等^[13]推导非饱和孔隙-裂隙膨胀土的外力功表达式，建立非饱和孔隙-裂隙膨胀土的弹塑性本构模型. Guo等^[14]在Borja等^[12]理论的基础上利用2个固相骨架应变和共轭的2个有效应力建立非饱和孔隙-裂隙土的双有效应力弹塑性模型. 基于嵌套思路，胡亚元^[15]建立考虑岩体（固相）材料变形的饱和孔隙-裂隙介质的本构框架，该模型选用的有效应力公式与Terzaghi公式存在较大差异.

本研究1）在小应变条件下，明确孔隙骨架体应变、裂隙骨架体应变、岩块材料体应变的物理含义；基于混合物理论中质量守恒方程，推导固相和流相的应变分解式，揭示流相与固相之间的耦合机理；为了更好地模拟岩体内部各组分之间的相互作用，考虑了更符合工程实际的岩体固相材料变形. 2）将固相骨架体应变、偏斜应变以及固流两相的材料体应变选作本构方程的应变状态变量，推导饱和孔隙-裂隙岩体的能量守恒方程，确定建立本构模型所需的共轭应力状态变量. 3）根据自由能势函数和热力学共轭力学量，建立饱和孔隙-隙岩体的线弹性本构方程，推导模型参数的取值方法.

如图1所示，饱和孔隙-裂隙岩体是由岩块固相S、裂隙流相F和孔隙流相P组成的混合物，固体和流体均可压缩. 由于岩体中孔隙和裂隙相互连通，在压力差作用下，孔隙和裂隙中部分流体从孔隙流入裂隙或从裂隙流入孔隙. 假定裂隙流相与孔隙流相间存在质量交换，固相与裂隙和孔隙中的流相不存在质量交换. 令$\alpha \in \{ {\text{S,F,P}}\} $为组分指代变量，${\varphi _\alpha }$为第$\alpha $组分的体积分数，表达式为${\varphi _\alpha } = {V_\alpha }/V$，其中${V_\alpha }$为第$\alpha $组特征单元体（representative elementary volume, REV）的体积，$V$为饱和孔隙-裂隙岩体混合物总体积. 令${\rho _\alpha }$为在饱和孔隙-裂隙岩体中第$\alpha $组分的平均密度，${\rho _{{\text{R}}\alpha }}$为第$\alpha $组分的基质密度（或称真实密度），满足${\rho _\alpha } = {\varphi _\alpha }{\rho _{{\text{R}}\alpha }}$，饱和孔隙-裂隙岩体混合物的总密度为$\rho = {\rho _{\text{S}}}+{\rho _{\text{F}}}+{\rho _{\text{P}}}$. 假定孔隙和裂隙被流体组分全部填充，根据体积分数定义，在饱和孔隙-裂隙岩体中有

定义固相、孔隙流相和裂隙流相的物质导数为d^α（·）/dt. 假定饱和孔隙-裂隙岩体中固相与流相间不存在质量交换，裂隙流相与孔隙流相间存在质量交换，则固相、裂隙流相以及孔隙流相的质量守恒方程分别表示为

式中：t为时间；$ {{\boldsymbol{v}}_{\mathrm{S}}} $，$ {{\boldsymbol{v}}_{\mathrm{F}}} $，$ {{\boldsymbol{v}}_{\mathrm{P}}} $分别为饱和孔隙-裂隙岩体固相、裂隙流相、孔隙流相的绝对速度；$ {\hat c_{\mathrm{F}}} $，$ {\hat c_{\mathrm{P}}} $为裂隙流相和孔隙流相间的质量交换率，满足$ {\hat c_{\mathrm{F}}}+{\hat c_{\mathrm{P}}} = 0 $. 把固相作为参考构形，则裂隙流相相对于固相的扩散速度为${{\boldsymbol{\tilde v}}_{\text{F}}} = {{\boldsymbol{v}}_{\text{F}}} - {{\boldsymbol{v}}_{\text{S}}}$^[12]，孔隙流相相对于固相扩散速度为${{\boldsymbol{\tilde v}}_{\text{P}}} = {{\boldsymbol{v}}_{\text{P}}} - {{\boldsymbol{v}}_{\text{S}}}$^[12]. 把${{\boldsymbol{\tilde v}}_{\text{F}}}$、${{\boldsymbol{\tilde v}}_{\text{P}}}$和${\rho _\alpha } = {\varphi _\alpha }{\rho _{{\text{R}}\alpha }}$代入式(2)~式(4)得到

在饱和孔隙-裂隙岩体中，裂隙比定义为裂隙体积与固相体积的比值，表达式为${e_{\text{C}}} = {\varphi _{\text{F}}}/{\varphi _{\text{S}}}$；孔隙比定义为孔隙体积与固相体积的比值，表达式为$ {e_{\text{H}}} = {\varphi _{\text{P}}}/{\varphi _{\text{S}}} $. 体应变以体积膨胀为正，孔隙骨架体变形率、裂隙骨架体变形率、岩块材料体变形率分别为

令$ {{\boldsymbol{\varepsilon }}_{\text{S}}} $为固相应变张量，$ {\varepsilon _{{\text{SV}}}} $为固相体应变，在小应变条件下有$ {\dot \varepsilon _{{\text{SV}}}}{\text{ = }}\nabla \cdot {{\boldsymbol{v}}_{\text{S}}} $^[15]，将${e_{\text{C}}} = {\varphi _{\text{F}}}/{\varphi _{\text{S}}}$、$ {e_{\text{H}}} = {\varphi _{\text{P}}}/{\varphi _{\text{S}}} $、式(8)~式(10)代入式(5)，将计算结果对$t$积分可以得到

由式(11)可知，饱和孔隙-裂隙岩体固相体应变可以分解为裂隙骨架体应变、孔隙骨架体应变与岩块材料体应变.

同理，在小应变条件下，裂隙流相体应变$ {\dot \varepsilon _{{\text{FV}}}}{\text{ = }}\nabla \cdot {{\boldsymbol{v}}_{\text{F}}} $^[15]，孔隙流相体应变$ {\dot \varepsilon _{{\text{PV}}}}{\text{ = }}\nabla \cdot {{\boldsymbol{v}}_{\text{P}}} $^[15]. 将裂隙流相材料体应变$ {\varepsilon _{{\text{RFV}}}} = \ln\; ({\rho _{R{\text{F0}}}}/{\rho _{{\mathrm{RF}}}}) $、孔隙流体材料体应变$ {\varepsilon _{{\text{RPV}}}} = \ln\; ({\rho _{R{\text{P0}}}}/{\rho _{{\mathrm{RP}}}}) $对t求导后，与$ {\rho _{\text{F}}} = {\varphi _{\text{F}}}{\rho _{{\text{RF}}}} $、$ {\rho _{\text{P}}} = {\varphi _{\text{P}}}{\rho _{{\text{RP}}}} $以及式(11)代入式(3)和式(4)，分别将计算结果对$t$积分，得到

式(12)、式(13)分别是裂隙流相和孔隙流相体应变的分解式. 由式(11)、式(12)和式(13)可以看出，裂隙骨架体应变和孔隙骨架体应变同时影响固相体应变和流相体应变. 固相、孔隙流相以及裂隙流相之间存在变形耦合效应，这种耦合作用通过裂隙骨架与孔隙骨架之间的变形实现.

根据混合物理论，令${\boldsymbol{ \sigma}} $为饱和孔隙-裂隙岩体混合物的柯西应力张量，${{\boldsymbol{\sigma }}_\alpha }$为第$\alpha $组分的柯西应力张量，有

混合物中应力以拉为正，压力以压为正. 总压力$ {p_{\mathrm{T}}} = - {\boldsymbol{\sigma }}:{{\boldsymbol{I}} / 3} $，岩块材料所受的真实压力$ {p_{\mathrm{S}}} $与$ {{\boldsymbol{\sigma }}_{\mathrm{S}}} $满足$ {\varphi _{\text{S}}}{p_{\mathrm{S}}} = - {{\boldsymbol{\sigma }}_{\mathrm{S}}}:{{\boldsymbol{I}} / 3} $，裂隙孔压$ {p_{\mathrm{F}}} $与$ {{\boldsymbol{\sigma }}_{\mathrm{F}}} $满足$ {{\boldsymbol{\sigma }}_{\text{F}}} = - {\varphi _{\text{F}}}{p_{\text{F}}}{\boldsymbol{I}} $，孔隙孔压$ {p_{\text{P}}} $与$ {{\boldsymbol{\sigma }}_{\mathrm{P}}} $满足$ {{\boldsymbol{\sigma }}_{\text{P}}} = - {\varphi _{\text{P}}}{p_{\text{P}}}{\boldsymbol{I}} $. 利用上述关系式和式(14)得到

令${\xi _\alpha }$为第$\alpha $组分的内能密度，由于不考虑温度场，可以忽略热流项和热源项，能量守恒方程只有机械功，$\displaystyle {\sum}_\alpha {{\rho _\alpha }} {\dot \xi _\alpha } = {\sum}_\alpha {{{\boldsymbol{\sigma }}_\alpha }:{\text{ }}} {{\boldsymbol{\dot \varepsilon }}_\alpha } $. 利用式(14)，$ {\dot {\boldsymbol{\varepsilon}} _{\alpha {\text{V}}}} = {{\mathrm{div}}} {{\boldsymbol{v}}_\alpha } $、$ {{\boldsymbol{\sigma }}_{\text{F}}} = - {\varphi _{\text{F}}}{p_{\text{F}}}{\boldsymbol{I}} $、$ {{\boldsymbol{\sigma }}_{\text{P}}} = - {\varphi _{\text{P}}}{p_{\text{P}}}{\boldsymbol{I}} $，${{\boldsymbol{\tilde v}}_{\text{F}}} = {{\boldsymbol{v}}_{\text{F}}} - {{\boldsymbol{v}}_{\text{S}}}$、${{\boldsymbol{\tilde v}}_{\text{P}}}= {{\boldsymbol{v}}_{\text{P}}} - {{\boldsymbol{v}}_{\text{S}}}$，能量守恒方程表示为

将式(16)等号右边第一项用总平均应力$ {\sigma _{\text{m}}} = {\boldsymbol{\sigma }}:{\boldsymbol{I}}/3 $和体应变$ {\varepsilon _{{\text{SV}}}} $，偏斜应力张量$ {\boldsymbol{S}} = {\boldsymbol{\sigma }} - {\sigma _{\text{m}}}{\boldsymbol{I}} $和偏斜应变张量$ {{\boldsymbol{e}}_{\text{S}}} = {{\boldsymbol{\varepsilon }}_{\text{S}}} - ({\varepsilon _{{\text{SV}}}}/{\text{3}}){\boldsymbol{I}} $表示；右边第二项和第三项分别用式(6)和式(7)替换. 再将$ {\varepsilon _{{\text{RPV}}}} = \ln \; ({\rho _{{\mathrm{R}}{\text{P0}}}}/{\rho _{{\mathrm{RP}}}}) $，$ {\varepsilon _{{\text{RFV}}}} = \ln \; ({\rho _{{\text{RF0}}}}/{\rho _{{\mathrm{RF}}}}) $对t求导后，联合${e_{\text{C}}} = {\varphi _{\text{F}}}/{\varphi _{\text{S}}}$、$ {e_{\text{H}}} = {\varphi _{\text{P}}}/{\varphi _{\text{S}}} $、式(8)~式(11)、$ {\sigma _{\text{m}}} = - ({\varphi _{\text{F}}}{p_{\text{F}}}+ {\varphi _{\text{P}}}{p_{\text{P}}}+{\varphi _{\text{S}}}{p_{\text{S}}}) $一起代入上述替换所得式，有

式中：$ {\text{ }}f \in \{ {\text{F,P}}\} $，$ {\tilde \sigma _{{\text{Cm}}}} = {\sigma _{\text{m}}}+{p_{\text{F}}} $、$ {\tilde \sigma _{{\text{Hm}}}} = {\sigma _{\text{m}}}+{p_{\text{P}}} $分别为裂隙平均有效应力和孔隙平均有效应力. 由式(17)可知，偏斜应力$ {\boldsymbol{S}} $和偏斜应变${{\boldsymbol{e}}_{\text{S}}}$为功共轭对，裂隙平均有效应力$ {\tilde \sigma _{{\text{Cm}}}} $与裂隙骨架体应变${\varepsilon _{{\text{CV}}}}$为功共轭对，孔隙平均有效应力${\tilde \sigma _{{\text{Hm}}}}$与孔隙骨架体应变${\varepsilon _{{\text{HV}}}}$为功共轭对，各组分真实压力${p_\alpha }$与各组分材料体应变${\varepsilon _{{\text{R}}\alpha {\text{V}}}}$为功共轭对.

令$ {\xi ^*} = \displaystyle\sum\nolimits_\alpha {{\rho _\alpha }{\xi _\alpha }} $，在小应变条件下,$ {\varphi _\alpha } \approx {\varphi _{\alpha {\text{0}}}} $,$ {\rho _\alpha } \approx {\rho _{\alpha 0}} $，略去高次项有$ {\dot \xi ^*} = {{\mathrm{d}}} \left(\displaystyle\sum\nolimits_\alpha {{\rho _{\alpha 0}}{\xi _\alpha }}\right) / {\text{d}}t = \displaystyle\sum\nolimits_\alpha {\rho _{\alpha 0}} ({\text{d}}{\xi _\alpha }/ {\text{d}}t) = \displaystyle\sum\nolimits_\alpha {{\rho _\alpha }({{\text{d}}^\alpha }{\xi _\alpha }/{\text{d}}t)} $，代入式(17)得到

根据热力学局部平衡假定，式(18)中的内能可以表示为$ {\xi ^*} = {\xi ^*}(\eta ,{{\boldsymbol{e}}_{\text{S}}},{\varepsilon _{{\text{CV}}}},{\varepsilon _{{\text{HV}}}},{\varepsilon _{{\text{R}}\alpha {\text{V}}}}) $，对它求全微分得到

其中$ \eta $为熵，对比式(18)和式(19)：

引入亥姆霍茨自由能$ {\psi ^*}(\theta ,{{\boldsymbol{e}}_{\text{S}}},{\varepsilon _{{\text{CV}}}},{\varepsilon _{{\text{HV}}}},{\varepsilon _{{\text{R}}\alpha {\text{V}}}}) $，有$ {\psi ^*}(\theta ,{{\boldsymbol{e}}_{\text{S}}},{\varepsilon _{{\text{CV}}}},{\varepsilon _{{\text{HV}}}},{\varepsilon _{{\text{R}}\alpha {\text{V}}}}) = {\xi ^*}(\eta ,{{\boldsymbol{e}}_{\text{S}}},{\varepsilon _{{\text{CV}}}},{\varepsilon _{{\text{HV}}}},{\varepsilon _{{\text{R}}\alpha {\text{V}}}}) - \theta \eta ，$ 对该自由能势函数求全微分后代入式(20)得到

令饱和孔隙-裂隙岩体混合物的初始平衡状态为$ (\theta ,{{\boldsymbol{e}}_{\text{S}}},{\varepsilon _{{\text{CV}}}},{\varepsilon _{{\text{HV}}}},{\varepsilon _{{\text{R}}\alpha {\text{V}}}}) = ({\theta _0},0,0,0,0) $，在受到微小的扰动后，将达到新的平衡态$ (\theta +{\theta _\Delta },{{\boldsymbol{e}}_{\text{S}}},{\varepsilon _{{\text{CV}}}},{\varepsilon _{{\text{HV}}}},{\varepsilon _{{\text{R}}\alpha {\text{V}}}}) $，其中$ {\theta _\Delta } $为温度增量. 根据势函数本构方程的一般性质，若自由能函数是状态变量的二次多项式函数，可获得线弹性本构关系. 在材料各向同性假设下，根据线性热应力理论，温度变化只会引起线弹性体的体积变化. 因此，可不考虑固相体应变和温度变化对固相偏斜应变的影响，并且在饱和孔隙-裂隙介质研究中，通常认为各组分材料只与其自身的变形和温度有关，而与其他组分的变形无关. 利用这一假定，岩块材料的变形与流体材料的变形相互解耦. 根据混合物理论互易定理，固相的裂隙骨架应变、孔隙骨架应变和岩块材料应变之间相互解耦. 也就是说，本构方程中的自变量之间相互独立. 假定整个过程温度不变，即$ {\theta _\Delta } = 0 $，则小应变线弹性条件下亥姆霍茨自由能取为

将(22)代入式(21)得到

式中：$ {G_{\text{b}}} $为饱和孔隙-裂隙岩体的剪切模量；$ {C_{{\text{HH}}}} $、$ {C_{{\text{CC}}}} $、$ {C_{\alpha \alpha }} $分别为孔隙骨架、裂隙骨架、岩块材料以及流体材料的体积柔度系数. 通过对文献[6]中式(4a)、式(4b)等号两边求迹，得到式(23b)、式(23c). 同样地，当不考虑温度影响时，对文献[9]中式(12a)、式(12b)等式两边求迹，得到与式(23b)、式(23c)相同的结果. 文献[6]和文献[9]中的固相材料本构方程与本研究当“$ {\text{ }}\alpha = {\text{S }} $”时的式(23d)不同，它们没有考虑固相真实压力与固相材料压力之间的关系. 将$ {p_{\text{T}}} = - {\sigma _{\text{m}}} $、$ {\tilde \sigma _{{\text{Cm}}}} = {\sigma _{\text{m}}}+{p_{\text{F}}} $、$ {\tilde \sigma _{{\text{Hm}}}} = {\sigma _{\text{m}}}+ {p_{\text{P}}} $、式(15)和式(23b)~式(23d)代入式(11)得到

式中：$ {\beta _{\text{p}}} = ({C_{\text{M}}} - {\varphi _{{\text{S0}}}}{C_{{\text{SS}}}})/{C_{\text{b}}} $，$ {\text{ }}{\beta _{\text{F}}} = 1 - {C_{\text{M}}}/{C_{\text{b}}} $，$ {C_{\text{b}}} = {C_{{\text{CC}}}}+{C_{{\text{HH}}}}+{C_{{\text{SS}}}} $，$ {C_{\text{M}}} = {C_{{\text{HH}}}}+{\varphi _{{\text{SP0}}}}{C_{{\text{SS}}}} $. 定义裂隙流体渗入量为${\zeta _{\text{F}}} = {\varphi _{{\text{F0}}}}({\varepsilon _{{\text{SV}}}} - {\varepsilon _{{\text{FV}}}})$，孔隙流体渗入量为${\zeta _{\text{P}}} = {\varphi _{{\text{P0}}}}({\varepsilon _{{\text{SV}}}} - {\varepsilon _{{\text{PV}}}})$. 利用式(11)~式(13)、式(15)、式(23b)~式(23d)、式(24)、$ {p_{\text{F}}} - {p_{\text{P}}} = {\gamma _{\text{w}}}{\hat c_{\text{F}}}/{\overline \alpha } {k_{\text{P}}}{\rho _{{\text{RF}}}} $^[16]，则${\zeta _{\text{F}}}$和${\zeta _{\text{P}}}$的表达式分别为

式中：$\overline \alpha $为形状系数，$ {k_{\text{F}}} $、$ {k_{\text{P}}} $分别为裂隙和孔隙流体的渗透系数，${\gamma _w}$为流体的重度，

式(24)、式(25)、式(26)分别为小应变条件下饱和孔隙-裂隙岩体固相、裂隙流相渗入量以及孔隙流相渗入量的本构方程，结合达西定律$ {\dot \zeta _{\text{F}}} = {k_{\text{F}}}{\nabla ^2}{p_{\text{F}}}/{\gamma _w} $、$ {\dot \zeta _{\text{P}}} = {k_{\text{P}}}{\nabla ^2}{p_{\text{P}}}/{\gamma _w} $，则饱和孔隙-裂隙岩体固结控制方程为

本研究的饱和孔隙-裂隙岩体水力耦合模型包括以下参数：体积分数${\varphi _{{\text{S0}}}}$、${\varphi _{{\text{F0}}}}$和${\varphi _{{\text{P0}}}}$；渗透系数${k_{\text{F}}}$、${k_{\text{P}}}$，形状系数$\overline \alpha $；弹性柔度系数${C_{{\text{SS}}}}$、${C_{{\text{HH}}}}$、${C_{{\text{CC}}}}$、${C_{{\text{R}}f}}$，剪切模量$ {G_{\text{b}}} $. ${\varphi _{{\text{F0}}}}$、${\varphi _{{\text{P0}}}}$由压汞试验得到，${\varphi _{{\text{S0}}}}$由体积分数之和为1的条件间接求得. 考虑到实际工程中无法保证孔隙、裂隙完全充满汞液体，因此结合核磁共振频谱分析、CT扫描测量方法^[17]减小测试过程造成的误差，采用敏感性分析量化饱和度造成的误差影响程度. 孔隙-裂隙岩体的渗透系数由野外水文地质试验、反演分析、公式估算及室内测试等方法获取^[18-19]. 形状系数$\overline \alpha $通过文献[4]中公式计算得到.

在实际工程中，孔隙-裂隙岩体的裂隙分布比较稀疏，特征单元体一般为1~2 m，甚至更大，因此很难进行相应的室内试验确定模型的固相弹性柔度系数. 为此，借鉴文献[6]确定模型参数的思想，视孔隙-裂隙岩体为组合体：1）将岩块固相材料与孔隙组成岩块基质，在岩体力学中称为完整岩块；2）将完整岩块视为整体与节理裂隙一起构成孔隙-裂隙岩体^[20]；3）分别测量岩块和节理裂隙的力学参数.

根据文献[15]中式(14)得到饱和孔隙-裂隙岩体中完整岩块的质量守恒方程为

在小应变条件下，$ \dot \varepsilon _{{\text{SV}}}^{\text{M}} = \nabla \cdot {\boldsymbol{v}}_{\text{S}}^{\text{M}} $为饱和孔隙-裂隙岩体中完整岩块的体变形率. 将${e_{\text{C}}} = {\varphi _{\text{F}}}/{\varphi _{\text{S}}}$、$ {e_{\text{H}}} = {\varphi _{\text{P}}}/{\varphi _{\text{S}}} $、式(8)~式(10)、$ \dot \varepsilon _{{\text{SV}}}^{\text{M}} = \nabla \cdot {\boldsymbol{v}}_{\text{S}}^{\text{M}} $、$ {\dot \varepsilon _{{\text{SV}}}}{\text{ = }}\nabla \cdot {{\boldsymbol{v}}_{\text{S}}} $代入式(30)，并将计算结果对$t$积分，得到

式(31)表明完整岩块的体应变可以分解为岩块材料的体应变和孔隙骨架体应变.

文献[7]的思想实验认为：当围压与裂隙孔压相同（$ {p_{\text{F}}} = - {\sigma _{\text{m}}} $）时，饱和孔隙-裂隙岩体中完整岩块的弹性系数与单独的完整岩块的弹性系数相同. 一方面，将$ {p_{\text{F}}} = - {\sigma _{\text{m}}} $、式(23b)~式(23d)代入式(31)得到饱和孔隙-裂隙岩体中完整岩块体应变为

另一方面，对饱和完整岩块试样单独进行力学试验，测定排水条件下完整岩块的体积模量为$K_{\text{b}}^{\text{M}}$和Biot系数$\alpha _{\text{P}}^{\text{M}}$. 完整岩块的体积应变表示为

式中：$\alpha _{\text{P}}^{\text{M}} = 1 - {C_{{\text{RS}}}}K_{\text{b}}^{\text{M}}$；${C_{{\text{RS}}}}$为岩块材料的柔度系数，通过等应力试验测定^[21]. 对比式(32)和式(33)，得到

结合$ K_{\text{b}}^{\text{M}} = {E_{\text{r}}}/[3(1 - 2{\mu _{\text{r}}})] $和式(34)进一步推出

完整岩块的弹性模量$ {E_{\text{r}}} $和泊松比$ {\mu _{\text{r}}} $由单轴压缩试验测出. 比较式(35)和文献[6]中孔隙骨架体积柔度系数式(14a)可以看出，文献[6]将完整岩块的变形视为孔隙骨架的变形，没有考虑岩块材料的压缩变形（$ {C_{{\text{RS}}}} = 0 $），而在岩体中岩块材料的变形通常不能忽略^[22]. 此外，在裂隙较为发育、岩体较为破碎的地质环境，$ {\varphi _{{\text{SP0}}}} $的数值与1相差较大，若采用文献[6]提出的孔隙骨架柔度系数计算，计算误差可能较大.

在荷载作用下，饱和孔隙-裂隙岩体体应变${\varepsilon _{{\text{SV}}}}$由节理裂隙变形引起的体应变${\varepsilon _{{\text{jV}}}}$和完整岩块体应变$ \varepsilon _{{\text{SV}}}^{\text{M}} $（包括岩块材料的体应变）共同组成. 根据式(11)和式(31)，有

对于结构面组数较多，随机定向分布的孔隙-裂隙岩体，通常可以认为是各向同性介质^[20]. 因此，裂隙节理变形引起的岩体体应变本构关系为

节理裂隙的法向刚度$ {k_{\text{n}}} $可与剪切刚度$ {k_{\text{τ }}} $由室内试验剪切仪测定；节理裂隙平均间距$d$通过野外水文地质试验测量得到^[8]. 将式(37)和${\varepsilon _{{\text{CV}}}} = {C_{{\text{CC}}}}{\sigma _{\text{m}}}$、${\varepsilon _{{\text{HV}}}} = {C_{{\text{HH}}}}{\sigma _{\text{m}}}$代入式(36)得到

对比式(38)和文献[6]的式(14b)可以发现，文献[6]的裂隙骨架体积柔度系数中缺少与孔隙骨架体积变形有关项$ {\varphi _{{\text{F0}}}}{C_{{\text{HH}}}}/{\varphi _{{\text{SP0}}}} $. 对比结果表明在节理法向刚度测量试验中，文献[6]将岩块视为刚体，忽略了完整岩块中孔隙变形，只考虑了裂隙的压缩，与实际情况不符.

令饱和孔隙-裂隙岩体中完整岩块的应变张量为${\boldsymbol{\varepsilon }}_{\text{S}}^{\text{M}}$，则偏斜应变张量$ {\boldsymbol{e}}_{\text{S}}^{\text{M}} = {\boldsymbol{\varepsilon }}_{\text{S}}^{\text{M}} - \varepsilon _{{\text{SV}}}^{\text{M}}{\boldsymbol{I}}/3 $，偏斜应力张量为$ {{\boldsymbol{S}}_{\text{M}}} = {{\boldsymbol{\sigma }}_{\text{M}}} - {\sigma _{{\text{Mm}}}}{\boldsymbol{I}} = ({\boldsymbol{\sigma }} - {\sigma _{\text{m}}}{\boldsymbol{I}})/{\varphi _{{\text{SP0}}}} $，其中$ {{\boldsymbol{\sigma }}_{\text{M}}} = {{\boldsymbol{\sigma }}_{{\text{SP}}}}/{\varphi _{{\text{SP0}}}} $^[15-16]为饱和孔隙-裂隙岩体中完整岩块应力张量. 与3.1节类似，有

完整岩块的剪切模量$ G_{\text{b}}^{\text{M}} = {E_{\text{r}}}/[2(1+{\mu _{\text{r}}})] $. 令饱和孔隙-裂隙岩体中节理裂隙变形引起的岩体应变张量为${{\boldsymbol{\varepsilon }}_{\text{j}}}$，偏斜应变张量$ {{\boldsymbol{e}}_{\text{j}}} = {{\boldsymbol{\varepsilon }}_{\text{j}}} - {\varepsilon _{{\text{jV}}}}{\boldsymbol{I}}/3 $，偏斜应力张量为$ {{\boldsymbol{S}}_{\text{j}}} = {{\boldsymbol{\tilde \sigma }}_{\text{C}}} - {\sigma _{{\text{Cm}}}}{\boldsymbol{I}} = {\boldsymbol{\sigma }} - {\sigma _{\text{m}}}{\boldsymbol{I}} $，则

节理裂隙的剪切模量$ {G_{{\text{CC}}}} = 1/{{\text{k}}_{\text{τ }}}d $. 为了简便，假定节理结构面泊松比${\mu _{\text{j}}}$和岩块泊松比${\mu _{\text{r}}}$相同，$ {G_{{\text{CC}}}} = 3(1 - 2{\mu _{\text{r}}})/[2(1+{\mu _{\text{r}}}){C_{{\text{CC}}}}] $. 由式(36)和${{\boldsymbol{\varepsilon }}_{\text{S}}} = {{\boldsymbol{\varepsilon }}_{\text{j}}}+ {\boldsymbol{\varepsilon }}_{\text{S}}^{\text{M}}$，得到

将$ {\boldsymbol{S}} = {\boldsymbol{\sigma }} - {\sigma _{\text{m}}}{\boldsymbol{I}} $、式(23a)、式(39)、式(40)代入式(41)，有

煤作为特殊的沉积岩，具有典型的双重孔隙结构^[23]. 根据立方定律建立裂隙渗透率-应力间的非线性关系，采用Anderson 01煤和Gilson 02煤的相关试验数据^[24]验证本研究构建的模型的有效性，试验过程保持孔隙和裂隙气压稳态. 模型拟合结果和试验数据如图2所示，拟合采用的参数如表1所示. 由立方定律^[23]可以得到

等式(9)两边对t求积分，得到

将式(45)、$ {\varphi _{\text{F}}} = {V_{\text{F}}}/V $、$ {\varphi _{{\text{F0}}}} = {V_{{\text{F0}}}}/{V_{\text{0}}} $代入式(44)得到

再将式(23b)、式(24)、式(46)、$ {\tilde \sigma _{{\text{Cm}}}} = {\sigma _{\text{m}}}+{p_{\text{F}}} $，$ {p_{\text{F}}} = {p_{\text{P}}} $代入式(43)，得到

从图2可以看出，由于压密作用，裂隙体积分数的降低导致渗透率逐渐减小，模型能够较好地拟合试验数据.

通过饱和孔隙-裂隙岩体偏斜应变本构方程式(23a)和固相体应变本构方程式(24)，推导出一维完全侧限条件下固相本构方程为

式中：${E_{\text{S}}} = (4{G_{\text{b}}}{C_{\text{b}}}+3)/3{C_{\text{b}}} $. 竖向应力平衡方程为

由式(48)、式(49)推导出饱和孔隙-裂隙岩体的固相竖向应力平衡方程为

在一维条件下，$ {\varepsilon _{{\text{S}}z}} = \partial {u_{{\text{S}}z}}/\partial z $，由式(28)、式(29)得到饱和孔隙-裂隙岩体一维固结方程为

采用厚度$h$的饱和孔隙-裂隙岩体进行一维固结特性分析，顶面($z = 0$)为自由排水边界，施加瞬时均布荷载${\sigma _0}$，底面($z = h$)为刚性不透水边界，竖向渗透系数保持恒定. 边界条件为

初始条件为

初始孔压$ {p_{{\text{P0}}}} $、$ {p_{{\text{F0}}}} $以及初始竖向应变$ {\varepsilon _{{\text{S}}z{\text{0}}}} $在一维条件下由式(25)、式(26)、式(48)结合初始条件式(54)计算得到，方程组为

将式(50)两边对$z$求积分，再对时间$t$求导后，代入式(51)和式(52)得到

用分离变量法解偏微分方程，得到$ t > 0 $时的解析解为

式中：${r_1}$、${r_2}$为方程${b_{\text{0}}}{r^2}+{b_{\text{1}}}r+{b_{\text{2}}} = 0$的解，

结合边界条件式(53c)、式(53d)、式(58)、式(59)，将等式(50)两边在区间$(0,h)$对$z$求两次积分后，得到任意时刻$t$饱和孔隙-裂隙岩体的固结沉降量为

为了直观地展示参数取值方法的理论差异带来的影响，反映岩体固相材料和孔隙骨架变形的影响，以一维固结算例进行比较分析. 饱和孔隙-裂隙岩体深$h = 20\;{\text{m}}$，上覆均布荷载${\sigma _0} = 4\;{\text{MPa}}$，边界条件同式(53). 综合文献[6]、文献[9]，选取完整岩块和节理裂隙基本力学参数如表2所示. 用本研究的参数取值方法和文献[6]中参数取值方法分别计算饱和孔隙-裂隙岩体的力学参数，计算结果如表3所示. 对比表3发现，根据本研究推导的参数取值方法计算的孔隙骨架体积柔度系数C_HH比文献[6]的小，裂隙骨架体积柔度系数C_CC比文献[6]的略大. 因此，与文献[6]参数确定方法相比，按本研究推导的参数取值方法计算的孔隙骨架的压缩性更小，裂隙骨架压缩性略大.

固结过程中基底以下8 m处，按表2、表3中模型参数计算所得的孔隙超孔压和裂隙超孔压的变化情况如图3所示. 可以看出，文献[6]方法计算的初始裂隙超孔压和初始孔隙超孔压分别为3.71 MPa、3.05 MPa，均高于本研究方法计算的3.62 MPa、2.86 MPa. 本研究方法计算得到的裂隙超孔压在均布荷载施加后201 s开始消散，此时孔隙超孔隙压不仅没有下降，反而上升，直到350 s后逐步消散. Khalili^[25]认为这是由于裂隙中超孔压消散，荷载从裂隙骨架转移到完整岩块上，体现了孔、裂隙和流体间的微观耦合作用. 相比而言，按照文献[6]方法计算得到的孔隙和裂隙超孔压开始消散时间更接近.

如图4所示为2种取值方法对饱和孔隙-裂隙岩体固结沉降的影响. 沉降过程经历3个阶段：缓慢增长、快速增长和沉降稳定，主要受孔隙和裂隙中超孔压消散和固相骨架变形的影响. 固结初期，由于岩体内超孔压较高，且沿深度均匀分布，相邻岩层间超孔压梯度较小，流体很难从排水层排除，沉降处于缓慢增长阶段；随着时间推移，靠近排水层的岩体超孔压率先消散，岩层间超孔压梯度逐渐增大，流体在压力梯度的作用下更易从顶部排水层排除，沉降进入快速增长阶段，直到沉降稳定，固结完成. 由于本研究方法考虑了岩块材料变形的影响，计算的初始沉降量为4.5 mm，大于文献[6]方法计算的初始沉降量(3.03 mm)，但总沉降量二者几乎相同.

在饱和孔隙-裂隙岩体中，裂隙与孔隙渗透系数的相对大小不仅影响孔隙和裂隙中超孔压的消散，还会影响岩体变形随固结时间变化的规律. 流体交换项$ \overline \alpha {k_{\text{P}}}({p_{\text{F}}} - {p_{\text{P}}})/{\gamma _w} $（或$ \overline \alpha {k_{\text{P}}}({p_{\text{P}}} - {p_{\text{F}}})/{\gamma _w} $）与孔隙渗透系数有关，因此保持孔隙渗透系数不变，将裂隙渗透系数减小至原来1/10和增大至原来10倍来调节渗透系数相对大小. 如图5所示为不同渗透系数比下的超孔压消散图. 可以看到，增大${k_{\text{F}}}/{k_{\text{P}}}$，提高了裂隙相对于孔隙的排水能力，裂隙中流体流速加快，缩短了从岩体顶部排出时间，加快了裂隙超孔压消散. 对孔隙超孔压而言，当${k_{\text{F}}}/{k_{\text{P}}}$相对较小时，固结初期出现积聚升高的现象，这是由于孔隙渗透能力和裂隙渗透能力相对接近，部分裂隙流体流入孔隙；随着${k_{\text{F}}}/{k_{\text{P}}}$增加，裂隙流体形成优先流，孔隙和裂隙间的超孔压力梯度加大，促使孔隙水外逸，孔隙超孔压消散加快，积聚值降低. 如图6所示为饱和孔隙-裂隙岩体固结沉降随$ {k_{\text{F}}}/{k_{\text{P}}} $的变化情况. 可以看到，随着$ {k_{\text{F}}}/{k_{\text{P}}} $增大，固结速率显著加快，固结完成时间更短.

形状系数与正交裂隙的组数和被切割岩块的大小有关，反映孔隙-裂隙中流体相互流动的特性，为此分析形状系数（$\overline \alpha $=1.5，15，150）对饱和孔隙-裂隙岩体固结性状的影响. 如图7所示为不同形状系数下超孔压消散图，当$\overline \alpha = 1.5$和$\overline \alpha = 15$时，裂隙超孔压消散曲线几乎重合，表明当形状系数较小时，增大形状系数对裂隙超孔压消散几乎没有影响，但对孔隙超孔压影响较大；随着形状系数增大，固结初期孔隙超孔压升高时间逐渐提前，平台期逐渐缩短. 这是由于在初始压力梯度的驱动下，裂隙流体逐渐流向孔隙中，孔隙超孔压上升；形状系数的增大使孔隙和裂隙间的连通性增强，缩短了孔隙和裂隙超孔压达到平衡的时间，在图中表现为实线和虚线开始重合的时间越来越早. 如图8所示，增大形状系数对岩体固结沉降的影响不明显.

（1）在小应变条件下，固流两相间的耦合作用通过裂隙骨架变形和孔隙骨架变形进行传递和协调；孔隙Terzaghi平均有效应力、裂隙Terzaghi平均有效应力、岩块材料真实压力、固相偏斜应力和流相材料真实压力是建立本构模型需要的应力状态变量.

（2）明确了孔隙、裂隙骨架变形的物理含义，强调了岩块材料的压缩性、岩块基质体积分数对孔隙弹性系数取值，以及孔隙骨架变形对裂隙弹性系数取值的影响，完善了参数的取值方法和理论依据. 一维固结算例表明，按本研究所提参数取值计算的超孔压消散规律以及初始沉降量与按文献[6]参数取值计算的结果差异较大，在岩体中，固相材料变形不可忽略.

（3）一维固结参数敏感性分析表明：1）增大渗透系数比会加快超孔压的消散，提高岩体的固结沉降速度；2）孔隙超孔压比裂隙超孔压对形状系数的变化更敏感，增大形状系数对岩体的固结沉降影响较小.

（4）本研究提出的模型为线性模型，后续计划针对饱和孔隙-裂隙岩体在极端条件下（如高温、高压）的非线性流固耦合行为展开研究.