文章快速检索     高级检索
  浙江大学学报(工学版)  2017, Vol. 51 Issue (2): 255-263  DOI:10.3785/j.issn.1008-973X.2017.02.005
0

引用本文 [复制中英文]

胡亚元. 非饱和多孔岩石的热力学本构理论[J]. 浙江大学学报(工学版), 2017, 51(2): 255-263.
dx.doi.org/10.3785/j.issn.1008-973X.2017.02.005
[复制中文]
HU Ya-yuan. Thermodynamics-based constitutive theory for unsaturated porous rock[J]. Journal of Zhejiang University(Engineering Science), 2017, 51(2): 255-263.
dx.doi.org/10.3785/j.issn.1008-973X.2017.02.005
[复制英文]

基金项目

国家自然科学基金资助项目(51178419)

作者简介

胡亚元(1968—)，男，副教授，从事地基处理和土的本构关系等研究.
orcid.org/0000-0002-5422-7679.
E-mail: huyayuan@zju.edu.cn

文章历史

收稿日期：2016-02-01
Contents              Abstract             Full text              Figures/Tables              PDF
非饱和多孔岩石的热力学本构理论
胡亚元     
浙江大学 滨海和城市岩土工程研究中心 浙江 杭州 310027
收稿日期：2016-02-01
基金项目：国家自然科学基金资助项目(51178419)
作者简介：胡亚元(1968—)，男，副教授，从事地基处理和土的本构关系等研究.
orcid.org/0000-0002-5422-7679.
E-mail: huyayuan@zju.edu.cn
摘要: 为了在本构模型中同时考虑岩石材料和多孔岩石的非线性和不可逆变形, 提出非饱和多孔岩石工程力学理论.应用体积分数概念和混合物理论, 凭借均匀化响应原理, 获得由孔隙应变张量、饱和度和各组分材料体应变5个状态变量表示的能量平衡方程.利用自由能确定的5个弹性方程, 加上体积分数之和等于1这个方程组成6个本构方程.根据这些方程可以求解非饱和岩石本构模型的全部6个未知变量(3个位移矢量和3个体积分数).根据不可逆热力学理论, 基于熵产公式提出用内变量表示的耗散率势函数, 获得能够反映黏性和塑性不可逆变形特性的耗散本构方程.结果表明, 自由能和耗散率2个势函数分别反映了岩石弹性和非弹性变形规律, 共同构成了非饱和多孔岩石的热力学本构理论框架.
关键词: 非饱和多孔岩石    平衡方程    熵产    Bishop有效应力    材料球应力    
Thermodynamics-based constitutive theory for unsaturated porous rock
HU Ya-yuan     
Research Center of Coastal and Urban Geotechnical Engineering, Zhejiang University, Hangzhou, 310058, China
Abstract: An engineering mechanics of unsaturated porous rock was proposed in order to consider the nonlinear and irreversible deformations of both rock material and porous rock. Firstly, with using the volume fraction concept and mixture theory, the equation of energy balance in terms of five state variables was obtained in virtue of the mixture homogenous response principle, which were called void strain tensor, degree of saturation and the material volume strain of each constituent. The five elastic equations which were determined by free energy and the one equation that the sum of all volume fractions was equal to unity compose six constitutive equations. Based on these equations, all six unknown variables (three displace vectors and three volume fractions) could be solved in the constitutive model of unsaturated rock. Secondly, according to the irreversible thermodynamics, the potential of dissipative rate in terms of internal variables was proposed on the basis of entropy production formula. Thus, the dissipative constitutive equations were derived to describe the irreversible deformation behaviors such as viscosity or plasticity. Results show that, the two potential functions of free energy and dissipative rate can depict the discipline of elastic and inelastic deformations, respectively, which both form the constitutive theoretical framework of unsaturated porous rock.
Key words: unsaturated porous rock    balance equation    entropy production    Bishop effective stress    material hydrostatic stress    

在岩土力学中, 把裂隙被液气2种流体充满的岩土称为非饱和多孔介质.随着大型水利水电工程和油气田开采工程的发展, 涉及裂隙岩石、液体和气体三相之间相互作用的工程问题越来越多, 促使近30年来非饱和多孔介质理论迅速发展[1-13].本文把固相颗粒变形称为固体材料变形.由于在多孔岩体中不能忽略岩石材料的变形, 因此无法把非饱和土理论完全应用于非饱和岩石领域[10-13].Buham等[11-15]建议采用与Skempton[9]公式类似的有效应力代替Bishop[16]有效应力来建立非饱和岩石力学理论, 但是, 类似Skempton公式的有效应力建立在线弹性理论基础之上, 它不适用于产生非线性或塑性变形的岩石材料[11, 13, 17].如何在考虑岩石材料变形的条件下建立非饱和岩石的一般本构理论体系, 到目前为止依然是一个有待解决的课题, 目前越来越多的岩土专家从热力学角度来研究非饱和岩石的本构理论[1][5][14-15].陈正汉等[14]建立了非饱和多孔岩土介质的新有效应力公式.赵成刚等[18]从组分连续条件出发推导了考虑固相材料变形的非饱和多孔介质机械功公式.Borja[16]研究了非饱和多孔介质细观和宏观力学特性之间的统计关系, 建立了考虑固体材料变形的非饱和多孔介质力学理论.但Borja[16]在推导本构方程时, 不合理地确定了体积分数本构方程, 致使把吸力所做的功全部视为由不可逆熵引起的耗散能, 结论显得欠合理.本文从传统混合物理论出发, 通过引入组分体积分数, 根据混合物理论的均匀化响应原理, 把非饱和多孔介质的能量方程表示成孔隙应变张量、饱和度和各组分材料体应变等广延量的微分表达式, 利用能量方程中的功共轭量建立自由能势函数, 以反映非饱和多孔介质的弹性特性；根据熵产公式结合内变量概念建立耗散率势函数, 以反映非饱和多孔介质的黏性和塑性等不可逆变形特性.建立了一个可以同时反映各组分可逆和不可逆变形的非饱和多孔介质本构理论框架.

1 建立自由能势函数 1.1 质量守恒

在符号上下标中, S表示非饱和岩土的固相, L表示液相和G表示气相.α∈{S, L, G}为组分指征变量.设第α组分的初始构型坐标为Xα, 时间t时的坐标为x, 则有

$ \mathit{\boldsymbol{x=}}{{\mathit{\boldsymbol{x}}}_{\alpha }}\left( {{\mathit{\boldsymbol{X}}}_{\alpha }},t \right)或{{\mathit{\boldsymbol{X}}}_{\alpha }}={{\mathit{\boldsymbol{X}}}_{\alpha }}\left( \mathit{\boldsymbol{x}},t \right). $ (1)

设φBα为第α组分的体积分数, 其定义为第α组分的体积占总体积之比.ρα、ρα分别为第α组分的密度和材料密度(也称真实密度), ρα=φBαρα.多孔岩石密度为$\rho = \sum\limits_\alpha {{\rho _\alpha }} $.令vα为第α组分的速度, 则非饱和多孔岩石的平均速度为

$ \mathit{\boldsymbol{v=}}\frac{1}{\rho }\sum\limits_{\alpha }{{{\rho }_{\alpha }}{{\mathit{\boldsymbol{v}}}_{\alpha }}}. $ (2)

第α组分相对于混合物平均速度的扩散速度用uα来表示, 其定义为

$ {{\mathit{\boldsymbol{u}}}_{\alpha }}={{\mathit{\boldsymbol{v}}}_{\alpha }}-\mathit{\boldsymbol{v}}. $ (3)

令Γα或Γα是定义在(x, t)上的场函数, 它关于第α组分运动的物质导数为[19]

$ \overset{\backslash }{\mathop{{{\mathit{\Gamma }}_{\alpha }}}}\,=\frac{\partial {{\mathit{\Gamma }}_{\alpha }}}{\partial t}+\rm{grad}{{\mathit{\Gamma }}_{\alpha }}\cdot {{\mathit{\boldsymbol{v}}}_{\alpha }}\ 或\ \overset{\backslash }{\mathop{{{\mathit{\Gamma }}^{\alpha }}}}\,=\frac{\partial {{\mathit{\Gamma }}^{\alpha }}}{\partial t}+\rm{grad}{{\mathit{\Gamma }}^{\alpha }}\cdot {{\mathit{\boldsymbol{v}}}_{\alpha }}. $ (4)

关于非饱和多孔岩石平均速度的物质导数为[15]

$ {{{\mathit{\dot{\Gamma }}}}_{\alpha }}\rm{=}\frac{\partial {{\mathit{\Gamma }}_{\alpha }}}{\partial t}+\rm{grad}{{\mathit{\Gamma }}_{\alpha }}\cdot \mathit{\boldsymbol{v}}\ 或\ {{{\mathit{\dot{\Gamma }}}}^{\alpha }}=\frac{\partial {{\mathit{\Gamma }}^{\alpha }}}{\partial t}+\rm{grad}{{\mathit{\Gamma }}^{\alpha }}\cdot \mathit{\boldsymbol{v}}. $ (5)

设φB=φBL+φBG为孔隙率, 饱和度Sr的定义为

$ {{S}_{\text{r}}}=\frac{{{\varphi }_{\text{BL}}}}{{{\varphi }_{\text{BL}}}+{{\varphi }_{\text{BG}}}}\times 100%. $ (6)

根据各组分的定义, 有

$ {{\varphi }_{\rm{BS}}}\left( \mathit{\boldsymbol{x}} \right)+{{\varphi }_{\rm{BL}}}\left( \mathit{\boldsymbol{x}} \right)+{{\varphi }_{\rm{BG}}}\left( \mathit{\boldsymbol{x}} \right)={{\varphi }_{\rm{BS}}}\left( \mathit{\boldsymbol{x}} \right)+{{\varphi }_{\rm{B}}}\left( \mathit{\boldsymbol{x}} \right)=1. $ (7)

假定非饱和多孔岩石各相之间不发生质量交换, 则各相的质量守恒条件为[19]

$ {{\overset{\backslash }{\mathop{\rho }}\,}_{\rm{S}}}+{{\rho }_{\beta }}\rm{div}\ {{\mathit{\boldsymbol{v}}}_{\rm{S}}}=0. $ (8)
$ {{\overset{\backslash }{\mathop{\rho }}\,}_{\beta }}+{{\rho }_{\beta }}\rm{div}\ {{\mathit{\boldsymbol{v}}}_{\beta }}=0,\ \ \ \ \beta \in \left\{ \rm{L},\rm{G} \right\}. $ (9)

式中：β∈{L, G}为液气2相流体的组分指征变量.根据式(6)-(7) 并利用式(2) 得

$ \dot{\rho }=+\rho \rm{div}\ \mathit{\boldsymbol{v=}}0. $ (10)

当Γα是第α组分的场函数时, Γ是Γα的质量加权平均值, 它的定义为

$ \mathit{\Gamma =}\frac{1}{\rho }\sum\limits_{\alpha }{{{\rho }_{\alpha }}{{\mathit{\Gamma }}_{\alpha }}}. $ (11)

根据式(2)、式(4)-(5) 和式(8)~(11) 可得[19]

$ \rho \mathit{\dot{\Gamma }=}\sum\limits_{\alpha }{\left[ {{\rho }_{\alpha }}\overset{\backslash }{\mathop{{{\mathit{\Gamma }}_{\alpha }}}}\,-\rm{div}\left( {{\rho }_{\alpha }}{{\mathit{\Gamma }}_{\alpha }}{{\mathit{\boldsymbol{u}}}_{\alpha }} \right) \right]}. $ (12)

把裂隙岩石现时位置作为混合物的参考构型, 则液气2相相对裂隙岩石现时位置的扩散速度为

$ {{\mathit{\boldsymbol{W}}}_{\rm{L}}}={{\mathit{\boldsymbol{v}}}_{\rm{L}}}-{{\mathit{\boldsymbol{v}}}_{\rm{S}}};{{\mathit{\boldsymbol{W}}}_{\rm{G}}}={{\mathit{\boldsymbol{v}}}_{\rm{G}}}-{{\mathit{\boldsymbol{v}}}_{\rm{S}}}. $ (13)

把ρα=φBaρα代入到式(8)-(9), 利用式(4) 和(13) 得

$ \frac{{{\varphi }_{\rm{BS}}}}{{{\rho }^{\rm{S}}}}\overset{\backslash }{\mathop{{{\rho }^{\rm{S}}}}}\,+\overset{\backslash }{\mathop{{{\varphi }_{\rm{BS}}}}}\,+{{\varphi }_{\rm{BS}}}\rm{div}\ {{\mathit{\boldsymbol{v}}}_{\rm{S}}}=0. $ (14)
$ \frac{{{\varphi _{{\rm{B}}\beta }}}}{{{\rho ^\beta }}}\mathop {{\rho ^\beta }}\limits^\backslash + \mathop {{\varphi _{{\rm{B}}\beta }}}\limits^\backslash + {\varphi _{{\rm{B}}\beta }}{\rm{div}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} + {\varphi _{{\rm{B}}\beta }}{\rm{div}}\;{\mathit{\boldsymbol{W}}_\beta } = 0,\beta \in \left\{ {{\rm{L}},{\rm{G}}} \right\}. $ (15)
1.2 自由能势函数建立过程

设σα为第α组分的柯西应力, σ为非饱和多孔岩石总柯西应力.根据混合物理论它的定义为[19]

$ \mathit{\boldsymbol{\sigma = }}\sum\limits_\alpha {{\mathit{\boldsymbol{\sigma }}_\alpha }} . $ (16)

令${{\mathit{\boldsymbol{\hat p}}}_a}$为第α组分的动量供应量, bα为外体力密度, 各组分的动量守恒方程可表示为[19]

$ {\rho _{\rm{S}}}\mathop {{\mathit{\boldsymbol{v}}_{\rm{S}}}}\limits^\backslash = {\rm{div}}\;{\mathit{\boldsymbol{\sigma }}_{\rm{S}}} + {\rho _{\rm{S}}}{\mathit{\boldsymbol{b}}_{\rm{S}}} + {{\mathit{\boldsymbol{\hat p}}}_{\rm{S}}}. $ (17)
$ {\rho _\beta }\mathop {{\mathit{\boldsymbol{v}}_\beta }}\limits^\backslash = {\rm{div}}\;{\mathit{\boldsymbol{\sigma }}_\beta } + {\rho _\beta }{\mathit{\boldsymbol{b}}_\beta } + {{\mathit{\boldsymbol{\hat p}}}_\beta },\;\;\;\;\beta \in \left\{ {{\rm{L}},{\rm{G}}} \right\}. $ (18)

利用式(17)-(18) 以及当Γα= $\mathit{\boldsymbol{v}}_a^\backslash $时的式(12) 得

$ \begin{array}{l} \rho \mathit{\boldsymbol{\dot v = }}\sum\limits_\alpha {{\rho _\alpha }\mathop {{\mathit{\boldsymbol{v}}_\alpha }}\limits^\backslash } - \sum\limits_\alpha {{\rm{div}}\left( {{\rho _\alpha }{\mathit{\boldsymbol{u}}_\alpha } \otimes {\mathit{\boldsymbol{u}}_\alpha }} \right)} = \\ \;\;\;\;\sum\limits_\alpha {\left[ {{\rm{div}}\left( {{\mathit{\boldsymbol{\sigma }}_\alpha } - {\rho _\alpha }{\mathit{\boldsymbol{u}}_\alpha } \otimes {\mathit{\boldsymbol{u}}_\alpha }} \right) + {\rho _\alpha }{\mathit{\boldsymbol{b}}_\alpha }} \right]} + \sum\limits_\alpha {{{\mathit{\boldsymbol{\hat p}}}_\alpha }} . \end{array} $ (19)

由于非饱和多孔岩石总动量守恒与${{\mathit{\boldsymbol{\hat p}}}_a}$无关, 故由式(19) 得[19]

$ \sum\limits_\alpha {{{\mathit{\boldsymbol{\hat p}}}_\alpha }} = 0,\;\;\;\;\alpha \in \left\{ {{\rm{S}},{\rm{L}},{\rm{G}}} \right\}. $ (20)

由各相的角动量守恒方程可得各相的柯西应力张量σα是对称张量.

设$\mathit{\boldsymbol{\mathscr{E}}}_a$, qα, rα和${{\hat \varepsilon }_\alpha }$分别是第α组分的内能、热流向量、外热供应量和对第α组分的能量供应量, 第α组分的能量守恒方程可表示为[19]

$ {\rho _{\rm{S}}}\mathop {{\mathscr{E}_{\rm{S}}}}\limits^\backslash = \;{\mathit{\boldsymbol{\sigma }}_{\rm{S}}}:{\rm{grad}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} - {\rm{div}}\;{\mathit{\boldsymbol{q}}_{\rm{S}}} + {\rho _{\rm{S}}}{r_{\rm{S}}} + {{\hat \varepsilon }_{\rm{S}}}. $ (21)
$ \begin{array}{l} {\rho _\beta }\mathop {{\mathscr{E}_\beta }}\limits^\backslash = \\ \;\;\;\;{\mathit{\boldsymbol{\sigma }}_\beta }:{\rm{grad}}\;{\mathit{\boldsymbol{v}}_\beta } - {\rm{div}}\;{\mathit{\boldsymbol{q}}_\beta } + {\rho _\beta }{r_\beta } + {{\hat \varepsilon }_\beta },\beta \in \left\{ {{\rm{S}},{\rm{L}},{\rm{G}}} \right\}. \end{array} $ (22)

式(21)-(22) 中的${{\hat \varepsilon }_\alpha }$满足[19]

$ \sum\limits_\alpha {\left( {{{\hat \varepsilon }_\alpha } + {\mathit{\boldsymbol{u}}_\alpha } \cdot {{\mathit{\boldsymbol{\hat p}}}_\alpha }} \right) \in 0,\alpha \in \left\{ {{\rm{S,L,G}}} \right\}} . $ (23)

流体(液体和气体)组分应力与孔压之间的关系为

$ {\mathit{\boldsymbol{\sigma }}_\beta } = - {\varphi _{{\rm{B}}\beta }}{P_\beta }{\bf{1}}. $ (24)

使用式(24)、(13) 和(16) 得

$ \begin{array}{l} \sum\limits_\alpha {{\mathit{\boldsymbol{\sigma }}_\alpha }:{\rm{grad}}\;{\mathit{\boldsymbol{v}}_\alpha }} = \\ \;\;\;\;\;\mathit{\boldsymbol{\sigma }}:{\rm{grad}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} - {P_{\rm{L}}}{\varphi _{{\rm{BL}}}}{\rm{div}}{\mathit{\boldsymbol{W}}_{\rm{L}}} - {P_{\rm{G}}}{\varphi _{{\rm{BG}}}}{\rm{div}}{\mathit{\boldsymbol{W}}_{\rm{G}}}. \end{array} $ (25)

定义裂隙中总流体的平均孔隙压力为

$ \bar P = P_{\rm{L}}^ * {S_{\rm{r}}} + P_{\rm{G}}^ * \left( {1 - {S_{\rm{r}}}} \right). $ (26)

式中：PL*和PG*为名义液体和气体孔隙压力, 它们满足

$ P_{\rm{L}}^ * = {P_{\rm{L}}},P_{\rm{G}}^ * = {P_{\rm{G}}}. $ (27)

定义Bishop有效应力${\mathit{\boldsymbol{\tilde \sigma }}}$为

$ \mathit{\boldsymbol{\bar \sigma }} = \mathit{\boldsymbol{\sigma }} + {P_{\rm{G}}}{\bf{1}} - {S_{\rm{r}}}\left( {{P_{\rm{G}}} - {P_{\rm{L}}}} \right){\bf{1}} = \mathit{\boldsymbol{\sigma }} + \bar P{\bf{1}}. $ (28)

定义吸力和有效吸力为

$ s = P_{\rm{G}}^ * - P_{\rm{L}}^ * ;\tilde s = {\varphi _{\rm{B}}}s. $ (29)

把式(15) 代入式(25), 再利用式(14)、(7) 推导出来的$\varphi _{{\text{BS}}}^\backslash + \varphi _{_{{\text{BL}}}}^\backslash + \varphi _{_{{\text{BG}}}}^\backslash - {\mathit{\boldsymbol{W}}_{{\rm{Lgrad}}}}{\varphi _{{\rm{BL}}}} - {\rm{ }}{\mathit{\boldsymbol{W}}_{{\rm{Ggrad}}}}{\varphi _{{\rm{BG}}}} = 0$和式(26)~(29) 得

$ \begin{array}{l} \sum\limits_\alpha {{\mathit{\boldsymbol{\sigma }}_\alpha }:{\rm{grad}}\;{\mathit{\boldsymbol{v}}_\alpha }} = \mathit{\boldsymbol{\tilde \sigma }}:{\rm{grad}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} + {\varphi _{{\rm{BS}}}}\bar P\frac{{\mathop {{\rho ^{\rm{S}}}}\limits^\backslash }}{{{\rho ^{\rm{S}}}}} + \\ \;\;\;\;\sum\limits_\beta {{\varphi _{{\rm{B}}\beta }}{P_\beta }\frac{{\mathop {{\rho ^\beta }}\limits^\backslash }}{{{\rho ^\beta }}} - \tilde s\mathop {{S_{\rm{r}}}}\limits^\backslash } + \sum\limits_\beta {{P_\beta }{\mathit{\boldsymbol{W}}_\beta } \cdot {\rm{grad}}{\varphi _{{\rm{B}}\beta }}} . \end{array} $ (30)

式(21)-(22) 相加并利用式(20)、(23)、(27) 和(30) 得

$ \begin{array}{l} \sum\limits_\alpha {{\rho _\alpha }\mathop {{\mathscr{E}_\alpha }}\limits^\backslash } = \mathit{\boldsymbol{\tilde \sigma }}:{\rm{grad}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} + {\varphi _{{\rm{BS}}}}\bar P\frac{{\mathop {{\rho ^{\rm{S}}}}\limits^\backslash }}{{{\rho ^{\rm{S}}}}} - \tilde s\mathop {{S_{\rm{r}}}}\limits^\backslash + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_\beta {{\varphi _{{\rm{B}}\beta }}{P_\beta }\frac{{\mathop {{\rho ^\beta }}\limits^\backslash }}{{{\rho ^\beta }}}} + \sum\limits_\alpha {\left( { - {\rm{div}}\;{\mathit{\boldsymbol{q}}_\alpha } + {\rho _\alpha }{r_\alpha }} \right)} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_\beta {{\mathit{\boldsymbol{W}}_\beta } \cdot \left( {{P_\beta }{\rm{grad}}{\varphi _{{\rm{B}}\beta }} - {{\mathit{\boldsymbol{\hat p}}}_\beta }} \right)} . \end{array} $ (31)

令PT=-σ:1/3和${{\tilde \sigma }_{\rm{m}}}$ =${\bf{\tilde{ \pmb{\mathsf{\tilde σ}} }}}$ :1/3, PT为非饱和多孔岩石的总球应力, ${{{\bf{\tilde{ \pmb{\mathsf{ σ}} }}}}_{\rm{m}}}$为Bishop有效球应力.根据式(28) 有${{{\bf{\tilde{ \pmb{\mathsf{ σ}} }}}}_{\rm{m}}}$ =P-PT.岩石材料的体应变为${\vartheta _{\text{s}}}$=ln (ρS/ρS0), 式中的ρS0为岩石的初始材料密度.注意到$\mathit{\boldsymbol{\sigma }}$, ${\mathit{\boldsymbol{\tilde \sigma }}}$和x以拉为正而PS和${\vartheta _{\text{s}}}$以压为正, 利用PT、${{\mathit{\boldsymbol{\tilde \sigma }}}_{\rm{m}}}$、${\vartheta _{\text{s}}}$的定义和它们的关系式, 式(31) 等式右边的前2项有

$ \begin{array}{l} \mathit{\boldsymbol{\tilde \sigma }}:{\rm{grad}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} + {\varphi _{{\rm{BS}}}}\bar P\frac{{\mathop {{\rho ^{\rm{S}}}}\limits^\backslash }}{{{\rho ^{\rm{S}}}}} = \\ \;\;\;\mathit{\boldsymbol{\tilde \sigma }}:\left( {{\rm{grad}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} + \frac{1}{3}\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash {\bf{1}}} \right) - \mathit{\boldsymbol{\tilde \sigma }}:\frac{1}{3}\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash {\bf{1}} + {\varphi _{{\rm{BS}}}}\bar P\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash = \\ \;\;\;\mathit{\boldsymbol{\tilde \sigma }}:\left( {{\rm{grad}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} + \frac{1}{3}\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash {\bf{1}}} \right) + {\varphi _{{\rm{BS}}}}\frac{{{P_{\rm{T}}} - {\varphi _{\rm{B}}}\bar P}}{{{\varphi _{{\rm{BS}}}}}}\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash . \end{array} $ (32)
$ {P_{\rm{S}}} = \frac{{{P_{\rm{T}}} - {\varphi _{\rm{B}}}\bar P}}{{{\varphi _{{\rm{BS}}}}}}. $ (33)

根据式(33) 和图 1(b)可知PS是固体材料(岩块)承受的真实压强, 本文把它称为固体的材料球应力, φBSPS称为固体材料的平均球应力.

图 1 单元体示意图 Fig. 1 Diagram of representative volume element

令${\vartheta _\beta }$ =ln(ρβ/ρβ0), β∈{L, G}, 其物理含义为流体材料的体应变, ρβ0为流体的初始材料密度.把式(32)-(33) 和${\vartheta _\beta }$代入到式(31) 得

$ \begin{array}{l} \sum\limits_\alpha {{\rho _\alpha }\mathop {{\mathscr{E}_\alpha }}\limits^\backslash } = \mathit{\boldsymbol{\tilde \sigma }}:\left( {{\rm{grad}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} + \frac{1}{3}\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash 1} \right) + {\varphi _{{\rm{BS}}}}{P_{\rm{S}}}\mathop {{S_{\rm{r}}}}\limits^\backslash + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_\beta {{\varphi _{{\rm{B}}\beta }}{P_\beta }\mathop {{\vartheta _\beta }}\limits^\backslash - \tilde s\mathop {{S_{\rm{r}}}}\limits^\backslash } + \sum\limits_\alpha {\left( { - {\rm{div}}\;{\mathit{\boldsymbol{q}}_\alpha } + {\rho _\alpha }{r_\alpha }} \right)} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_\beta {{\mathit{\boldsymbol{W}}_\beta } \cdot \left( {{P_\beta }{\rm{grad}}{\varphi _{{\rm{B}}\beta }} - {{\mathit{\boldsymbol{\hat p}}}_\beta }} \right)} . \end{array} $ (34)

为了在内能平衡方程的基础上引入自由能势函数, 需要把它们表示成相互独立状态变量的函数, 然而在式(34) 中, ${\tilde s}$和PS与PL和PG之间存在关系式式(27), 因此${\tilde s}$和PS与PL和PG之间不满足独立条件.为了使式(34) 中的状态变量相互独立, 需要在${\tilde s}$和PS与PL和PG之间引入2个拉格朗日乘子λL和λG分别作用式(27) 中的2个等式, 令

$ \begin{array}{l} {0_{P_{\rm{L}}^ * - {P_{\rm{L}}}}} \equiv \\ \;\;\;\;0 = P_{\rm{L}}^ * - {P_{\rm{L}}},{0_{P_{\rm{G}}^ * - {P_{\rm{G}}}}} \equiv 0 = P_{\rm{G}}^ * - {P_{\rm{G}}}. \end{array} $ (35)

由式(35) 把拉格朗日乘子算法写成微分形式得

$ {0_{P_{\rm{L}}^ * - {P_{\rm{L}}}}}\mathop {{\lambda _{\rm{L}}}}\limits^\backslash \equiv 0\;aaa\;{0_{P_{\rm{G}}^ * - {P_{\rm{G}}}}}\mathop {{\lambda _{\rm{G}}}}\limits^\backslash \equiv 0. $ (36)

把式(35) 和(36) 代入式(34) 得

$ \begin{array}{l} \sum\limits_\alpha {{\rho _\alpha }\mathop {{\mathscr{E}_\alpha }}\limits^\backslash } = \mathit{\boldsymbol{\tilde \sigma }}:\left( {{\rm{grad}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} + \frac{1}{3}\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash {\bf{1}}} \right) + {\varphi _{{\rm{BS}}}}{P_{\rm{S}}}\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash + \\ \;\;\;\sum\limits_\beta {{\varphi _{{\rm{B}}\beta }}{P_\beta }\mathop {{\vartheta _\beta }}\limits^\backslash - \tilde s\mathop {{S_{\rm{r}}}}\limits^\backslash } + {0_{P_{\rm{L}}^* - {P_{\rm{L}}}}}\mathop {{\lambda _{\rm{L}}}}\limits^\backslash + {0_{P_{\rm{G}}^* - {P_{\rm{G}}}}}\mathop {{\lambda _{\rm{G}}}}\limits^\backslash + \\ \;\;\;\sum\limits_\alpha {\left( { - {\rm{div}}\;{\mathit{\boldsymbol{q}}_\alpha } + {\rho _\alpha }{r_\alpha }} \right)} + \sum\limits_\beta {{\mathit{\boldsymbol{W}}_\beta } \cdot \left( {{P_\beta }{\rm{grad}}{\varphi _{{\rm{B}}\beta }} - {{\mathit{\boldsymbol{\hat p}}}_\beta }} \right)} . \end{array} $ (37)

根据拉格朗日乘子算法理论可知, 式(37) 中的${\tilde s}$和PS与PL和PG之间不再受式(27) 的约束, 因此在式(37) 中${\tilde s}$、PS、PL和PG之间是相互独立的.

利用混合物的均匀化响应原理来说明能量平衡方程表示成式(37) 的优点.混合物均匀化响应原理的内容为：当多相多孔混合物承受外荷载时, 若单元体中每一点的真实应变增量(或速率)相等, 则该多相介质单元体等效于单相均匀单元体, 即单元体内每一点处的真实应力增量(或加荷速率)也相等；反之也然.由于岩石裂隙中的流体只能承受孔压, 因此, 混合物均匀化响应原理在非饱和多孔岩石中主要应用于研究各组分球应力和体应变的相互作用, 如Geertsma[20]、Skempton[9]和陈正汉[14]应用这一原理建立了饱和以及非饱和岩土的有效应力公式.

假定混合物均匀化响应原理适用于由任意n相组成的多孔介质(当然也适用于三相非饱和多孔岩石), 那么如图 1所示, 根据它可以获得如下3个结论：1) 当单元体内任意一点的应力加荷速率均匀即$P_{_{\text{T}}}^\backslash = P_{_{\text{L}}}^\backslash = P_{\text{G}}$时, 由式(27) 和(26) 可知$\mathop {\bar P}\limits^\backslash = P_{_{\text{L}}}^\backslash = P_{\text{G}}^\backslash = P_{_{\text{T}}}^\backslash $, 再根据由式(28) 推导出的${{\mathit{\boldsymbol{\tilde \sigma }}}_{\rm{m}}} = \bar{P}-{{P}_{\text{T}}}$可得$\mathop {{{{\bf{ \pmb{\mathsf{\tilde σ}} }}}_{\rm{m}}}}\limits^\backslash = \mathop {\bar P}\limits^\backslash - P_{\rm{T}}^\backslash \ $ =0；同时根据式(29) 可知$\mathop {\tilde s}\limits^\backslash $ =0. 2) 当单元体内任意一点的应力加荷速率均匀即$P_{_{\text{T}}}^\backslash = P_{_{\text{L}}}^\backslash = P_{\text{G}}^\backslash$时, 由式(27) 和(33) 可得$P_{\text{s}}^\backslash = P_{_{\text{L}}}^\backslash = P_{\text{G}}^\backslash = P_{_{\text{T}}}^\backslash $；由于固体的材料变形不可忽略, 故在$P_{\text{s}}^\backslash $作用下将产生体应变速率$\vartheta _{_{\text{s}}}^\backslash $.根据混合物均匀化响应原理此时单元体内任意一点的体应变增速相等, 故孔隙中的液体、气体和整个单元体均产生体应变速率$\vartheta _{_{\text{s}}}^\backslash $.注意到xS以拉为正而${\vartheta _{\text{s}}}$以压为正, 本段落上文1) 和2) 说明在$\mathit{\boldsymbol{\tilde \sigma }}_{\rm{m}}^\backslash $ =0和$\mathop {\tilde s}\limits^\backslash $ =0但P≠0时式(31) 中的多孔介质体应变grad vS:1= $ - \vartheta _{_{\text{s}}}^\backslash \;$ ≠0.这就表明grad vS除了受和的影响外还必然受的影响.如果选择${\bf{\tilde{ \pmb{\mathsf{ σ}} }}}$和P分别作为grad vS和$\vartheta _{_{\text{s}}}^\backslash $的热力学共轭量, 则grad vS和$\vartheta _{_{\text{s}}}^\backslash $之间必然存在相互耦合作用, 这就增加了建立非饱和多孔岩石本构关系的困难.3) 如果选用${\rm{grad}}\;\;{\mathit{\boldsymbol{v}}_{\rm{S}}} + \frac{1}{3}\vartheta _{_{\rm{s}}}^\backslash \boldsymbol{1}$作为与${\mathit{\boldsymbol{\tilde \sigma }}}$对偶的热力学应变共轭量, 那么根据本段落分析可得当$\mathop {\mathit{\boldsymbol{\tilde \sigma }}}\limits^\backslash $ =0和$\mathop {\tilde s}\limits^\backslash \;$ =0时必有${\rm{grad}}\;\;{\mathit{\boldsymbol{v}}_{\rm{S}}} + \frac{1}{3}\vartheta _{_{\rm{s}}}^\backslash \boldsymbol{1} \equiv 0$.这就表明如果混合物的均匀化响应原理成立, 那么在弹性条件下, ${\rm{grad}}\;\;{\mathit{\boldsymbol{v}}_{\rm{S}}} + \frac{1}{3}\vartheta _{_{\rm{s}}}^\backslash \boldsymbol{1}$的值只与${\mathit{\boldsymbol{\tilde \sigma }}}$和${\tilde s}$有关而与固体的材料应力PS无关；根据热力学对称性, 式(37) 中的岩石材料体应变${\vartheta _{\text{s}}}$亦只跟与${\vartheta _{\text{s}}}$对偶的应力PS有关, 这无疑能大大简化考虑岩石材料变形时非饱和岩石本构关系的建模工作.故本文把式(31) 转换成式(37) 的形式进行自由能势函数研究.由于混合物的均匀化响应原理还有待试验的进一步验证, 因此本文只把混合物均匀化响应原理视为实际混合物变形特性的一种理想化近似.为了更普遍地反映实际混合物的变形特性, 下文的理论推导如不特别说明假定混合物均匀化响应原理可能不成立.

令${\mathit{\boldsymbol{F}}_{\rm{S}}} = \frac{{\partial {\mathit{\boldsymbol{x}}_{\rm{S}}}}}{{\partial {\mathit{\boldsymbol{X}}_{\rm{S}}}}},{J_{\rm{S}}} = {\rm{det}}[\frac{{\partial {\mathit{\boldsymbol{x}}_{\rm{S}}}}}{{\partial {\mathit{\boldsymbol{X}}_{\rm{S}}}}}],{J^{\rm{S}}} = \frac{{{\rho ^{{\rm{S}}0}}}}{{{\rho ^{\rm{S}}}}},{\mathit{\boldsymbol{C}}_{\rm{S}}} = \mathit{\boldsymbol{F}}_{\rm{S}}^{\rm{T}} \cdot {\mathit{\boldsymbol{F}}_{\rm{S}}},$${\mathit{\boldsymbol{E}}_{\rm{S}}} = \frac{1}{2}(\mathit{\boldsymbol{F}}_{\rm{S}}^{\rm{T}} \cdot \mathit{\boldsymbol{F}} - {\bf{1}}),$${{\mathit{\boldsymbol{\tilde T}}}_{\rm{S}}} = \mathit{\boldsymbol{F}}_{\rm{S}}^{ - 1} \cdot \mathit{\boldsymbol{\tilde \sigma }} \cdot \mathit{\boldsymbol{F}}_{\rm{S}}^{{\rm{ - T}}},$${\mathit{\boldsymbol{F}}_{\rm{H}}} = {({J^{\rm{S}}})^{ - 1/3}}{\mathit{\boldsymbol{F}}_{\rm{S}}},,$${\mathit{\boldsymbol{C}}_{\rm{H}}} = \mathit{\boldsymbol{F}}_{\rm{H}}^{\rm{T}} \cdot {\mathit{\boldsymbol{F}}_{\rm{H}}},$${\mathit{\boldsymbol{E}}_{\rm{H}}} = \frac{1}{2}(\mathit{\boldsymbol{F}}_{\rm{H}}^{\rm{T}} \cdot {\mathit{\boldsymbol{F}}_{\rm{H}}} - {\bf{1}}),$${{\mathit{\boldsymbol{\tilde T}}}_{\rm{H}}} = \mathit{\boldsymbol{F}}_H^{ - 1} \cdot \mathit{\boldsymbol{\tilde \sigma }} \cdot \mathit{\boldsymbol{F}}_{\rm{H}}^{ - {\rm{T}}},$$\mathit{\boldsymbol{E}}_{\rm{H}}^\backslash = \mathit{\boldsymbol{F}}_{\rm{H}}^{\rm{T}} \cdot ({\rm{grad}}\;\mathit{\boldsymbol{x}}_{\rm{S}}^\backslash + \frac{1}{3}\vartheta _{\rm{s}}^\backslash {\bf{1}}) \cdot {\mathit{\boldsymbol{F}}_{\rm{H}}}.$式中：FS为固相(多孔岩石)的变形梯度, JS为FS的行列式值, 也是固相的体积比, CS为固相(多孔岩石)的右Cauchy-Green张量, ES为固相(多孔岩石)Green应变, ${{\mathit{\boldsymbol{\tilde{T}}}}_{\rm{s}}}$为固相(多孔岩石)的Piola-Kirchhoff有效应力；JS是固相材料的体积比；FH为孔隙引起的变形梯度, CH为FH的右Cauchy-Green张量, EH为孔隙Green应变, ${{\mathit{\boldsymbol{\tilde{T}}}}_{\rm{H}}}$为孔隙的Piola-Kirchhoff有效应力.把它们代入式(37) 后式(37) 变为

$ \begin{array}{*{20}{l}} {\sum\limits_\alpha {{\rho _\alpha }\mathop {{\mathscr{E}_\alpha }}\limits^\backslash } = {{\mathit{\boldsymbol{\tilde T}}}_{\rm{H}}}:\mathop {{\mathit{\boldsymbol{E}}_{\rm{H}}}}\limits^\backslash + {\varphi _{{\rm{BS}}}}{P_{\rm{S}}}\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash - \tilde s\mathop {{S_{\rm{r}}}}\limits^\backslash + \sum\limits_\beta {{\varphi _{{\rm{B}}\beta }}{P_\beta }\mathop {{\vartheta _\beta }}\limits^\backslash } + }\\ {\;\;\;{0_{P_{\rm{L}}^* - {P_{\rm{L}}}}}\mathop {{\lambda _{\rm{L}}}}\limits^\backslash + {0_{P_{\rm{G}}^* - {P_{\rm{G}}}}}\mathop {{\lambda _{\rm{G}}}}\limits^\backslash + \sum\limits_\alpha {\left( { - {\rm{div}}\;{\mathit{\boldsymbol{q}}_\alpha } + {\rho _\alpha }{r_\alpha }} \right)} + }\\ {\;\;\;\sum\limits_\beta {{\mathit{\boldsymbol{W}}_\beta } \cdot \left( {{P_\beta }{\rm{grad}}{\varphi _{{\rm{B}}\beta }} - {{\mathit{\boldsymbol{\hat p}}}_\beta }} \right)} = }\\ {\;\;\;{\rho _{\rm{S}}}\left( {\frac{{{{\mathit{\boldsymbol{\tilde T}}}_{\rm{H}}}}}{{{\rho _{\rm{S}}}}}:\mathop {{\mathit{\boldsymbol{E}}_{\rm{H}}}}\limits^\backslash + \frac{{{P_{\rm{S}}}}}{{{\rho ^{\rm{S}}}}}\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash - \frac{{\tilde s}}{{{\rho _{\rm{S}}}}}\mathop {{S_{\rm{r}}}}\limits^\backslash } \right) + \sum\limits_\beta {{\rho _\beta }\frac{{{P_\beta }}}{{{\rho ^\beta }}}\mathop {{\vartheta _\beta }}\limits^\backslash } + }\\ {\;\;\;{0_{P_{\rm{L}}^* - {P_{\rm{L}}}}}\mathop {{\lambda _{\rm{L}}}}\limits^\backslash + {0_{P_{\rm{G}}^* - {P_{\rm{G}}}}}\mathop {{\lambda _{\rm{G}}}}\limits^\backslash + \sum\limits_\alpha {\left( { - {\rm{div}}\;{\mathit{\boldsymbol{q}}_\alpha } + {\rho _\alpha }{r_\alpha }} \right)} + }\\ {\;\;\;\sum\limits_\beta {{\mathit{\boldsymbol{W}}_\beta } \cdot \left( {{P_\beta }{\rm{grad}}{\varphi _{{\rm{B}}\beta }} - {{\mathit{\boldsymbol{\hat p}}}_\beta }} \right)} .} \end{array} $ (38)

在以往的饱和多孔介质和非饱和多孔介质中, 通常假定孔隙流体材料的本构关系与单相介质单独存在时的本构关系一致, 即流体的材料应变只跟该流体的材料应力有关.本文也采用这一简化假定.从式(38) 的右边也可以看出非饱和多孔岩石的内能可以分解成多孔固体$\mathscr{E}_{\text{S}}$、液体$\mathscr{E}_{\text{L}}$和气体$\mathscr{E}_{\text{G}}$三者之和, 它们的表达式为：$\mathscr{E}_{\text{S}}$(ηS, EH, Sr, ${{\vartheta }_{\text{S}}}$, EHp, Srp, $\vartheta _{\text{S}}^{p}$, λL, λG)、$\mathscr{E}_{\beta }$ (ηβ, ${{\vartheta }_{\beta }}$, $\vartheta _{\beta }^{p}$), β∈{L, G}, 式中ηα为各组分的熵, EHp, Srp, $\vartheta _{\text{S}}^{p}$和$\vartheta _{\beta }^{p}$是固体和流体相的内应变.把它们代入到式(38) 得

$ \begin{array}{l} \sum\limits_\alpha {{\rho _\alpha }\mathop {{\mathscr{E}_\alpha }}\limits^\backslash } = {\rho _{\rm{S}}}\left( {\frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial {\eta _{\rm{S}}}}}\mathop {{\eta _{\rm{S}}}}\limits^\backslash + \frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial {\mathit{\boldsymbol{E}}_{\rm{H}}}}}:\mathop {{\mathit{\boldsymbol{E}}_{\rm{H}}}}\limits^\backslash + \frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial {\vartheta _{\rm{S}}}}}\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash + } \right.\\ \left. {\frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial {S_{\rm{r}}}}}\mathop {{S_{\rm{r}}}}\limits^\backslash + \frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial \mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash } \right) + {\rho _{\rm{S}}}\left( {\frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial \vartheta _{\rm{S}}^{\rm{p}}}}\mathop {\vartheta _{\rm{S}}^{\rm{p}}}\limits^\backslash + \frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial S_{\rm{r}}^{\rm{p}}}}:\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash } \right) + \\ \sum\limits_\beta {{\rho _\beta }\left( {\frac{{\partial {\mathscr{E}_\beta }}}{{\partial {\eta _\beta }}}\mathop {{\eta _\beta }}\limits^\backslash + \frac{{\partial {\mathscr{E}_\beta }}}{{\partial {\vartheta _\beta }}}\mathop {{\vartheta _\beta }}\limits^\backslash + \frac{{\partial {\mathscr{E}_\beta }}}{{\partial \vartheta _\beta ^{\rm{p}}}}\mathop {\vartheta _\beta ^{\rm{p}}}\limits^\backslash } \right)} + \\ {\rho _{\rm{S}}}\left( {\frac{{\partial {\mathscr{E}_\beta }}}{{\partial {\lambda _{\rm{L}}}}}\mathop {{\lambda _{\rm{L}}}}\limits^\backslash + \frac{{\partial {\mathscr{E}_\beta }}}{{\partial {\lambda _{\rm{G}}}}}\mathop {{\lambda _{\rm{G}}}}\limits^\backslash } \right). \end{array} $ (39)

根据热力学局部平衡假定, 可得

$ \theta = \frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial {\eta _{\rm{S}}}}} = \frac{{\partial {\mathscr{E}_\beta }}}{{\partial {\eta _\beta }}},{{\mathit{\boldsymbol{\tilde T}}}_{\rm{H}}} = {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial {\mathit{\boldsymbol{E}}_{\rm{H}}}}},\tilde s = - {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial {S_{\rm{r}}}}}. $ (40)
$ {P_{\rm{S}}} = {\rho ^{\rm{S}}}\frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial {\vartheta _{\rm{S}}}}},{P_\beta } = {\rho ^\beta }\frac{{\partial {\mathscr{E}_\beta }}}{{\partial {\vartheta _\beta }}}. $ (41)
$ {0_{P_{\rm{L}}^* - {P_{\rm{L}}}}} = {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial {\lambda _{\rm{L}}}}},{0_{P_{\rm{G}}^* - {P_{\rm{G}}}}} = {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial {\lambda _{\rm{G}}}}}. $ (42)

式中:θ为饱和岩石的温度.令

$ \mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}} = - {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial \mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}},{{\tilde s}^{\rm{p}}} = - {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial S_{\rm{r}}^p}}. $ (43)
$ P_{\rm{S}}^{\rm{p}} = - {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{\rm{S}}}}}{{\partial \vartheta _{\rm{S}}^p}},P_\beta ^{\rm{p}} = - {\rho _\beta }\frac{{\partial {\mathscr{E}_\beta }}}{{\partial \vartheta _\beta ^{\rm{p}}}}. $ (44)

式(43) 和(44) 中的$\mathit{\boldsymbol{\tilde{T}}}_{\rm{H}}^{\rm{P}}$, ${{{\tilde{s}}}^{\text{P}}}$, PSp和Pβp是相应的内应力.式(40) 和(41) 表明孔隙的Piola-Kirchhoff有效应力${{{\mathrm{\tilde{T}}}}_{\text{H}}}$与孔隙Green应变EH共轭, 有效吸力${\tilde{s}}$与饱和度Sr共轭, 由于固体和液体的材料密度是真实密度, 故固体材料球应力PS与固体材料体应变${{\vartheta }_{\text{S}}}$共轭, 流体材料的孔压Pβ与流体材料体应变${{\vartheta }_{\beta }}$共轭[21]. ${\mathit{\boldsymbol{D}}_{\rm{S}}}{\rm{ = }}\frac{1}{2}{\rm{[grad }}\;\;{\mathit{\boldsymbol{v}}_{\rm{S}}}{\rm{ + (grad }}\;\;{\mathit{\boldsymbol{v}}_{\rm{S}}}{{\rm{)}}^{\rm{T}}}{\rm{]}}$, DS为多孔岩石的应变率.令${\mathit{\boldsymbol{D}}_{\rm{H}}} = {\mathit{\boldsymbol{D}}_{\rm{S}}} + \frac{1}{3}\vartheta _{\rm{s}}^{\rm{\backslash }}{\kern 1pt} {\bf{1}}$, 故有

$ \begin{array}{l} {\rm{tr}}{\mathit{\boldsymbol{D}}_{\rm{S}}} = \\ \;\;{\rm{tr}}\;{\mathit{\boldsymbol{D}}_{\rm{H}}} - \mathop {{\vartheta _{\rm{S}}}}\limits^\backslash \;和\;{\mathit{\boldsymbol{D}}_{\rm{S}}} - \frac{{{\rm{tr}{\mathit{\boldsymbol{D}}_{\rm{S}}}}}}{{\bf{3}}} = {\mathit{\boldsymbol{D}}_{\rm{H}}} - \frac{{{\rm{tr}{\mathit{\boldsymbol{D}}_{\rm{H}}}}}}{{\bf{3}}}{\bf{1}}. \end{array} $ (45)

由于${\rm{tr}}\;{\mathit{\boldsymbol{D}}_{\rm{S}}}{\rm{ = - }}\rho _{\rm{S}}^{\rm{\backslash }}/{\rho _{\rm{S}}} = - (\varphi _{_{{\rm{BS}}}}^\backslash /{\varphi _{{\rm{BS}}}}) - \mathop {{\vartheta _{\rm{S}}}}\limits^\backslash = \left[{\mathop e\limits^\backslash /\left( {1 + e} \right)} \right] - \mathop {{\vartheta _{\rm{S}}}}\limits^\backslash $, 其中e为多孔岩石的孔隙比.由式(45) 的第一式可知${\rm{tr}}\;\;{\mathit{\boldsymbol{D}_{\rm{H}}}} = - \varphi _{_{{\rm{BS}}}}^\backslash /{\varphi _{{\rm{BS}}}} = \mathop e\limits^\backslash /\left( {1 + e} \right)$, 本文把DH称为孔隙应变率.根据${{\mathit{\boldsymbol{\tilde T}}}_{\rm{H}}}$、EH、DS和DH的定义有

$ \begin{array}{l} {{\mathit{\boldsymbol{\tilde T}}}_{\rm{H}}}:\mathop {{\mathit{\boldsymbol{E}}_{\rm{H}}}}\limits^\backslash = \left( {\mathit{\boldsymbol{F}}_{\rm{H}}^{ - 1} \cdot \mathit{\boldsymbol{\tilde \sigma }} \cdot \mathit{\boldsymbol{F}}_{\rm{H}}^{ - {\rm{T}}}} \right):\mathop {{\mathit{\boldsymbol{E}}_{\rm{H}}}}\limits^\backslash = \\ \;\;\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{\tilde \sigma }}:\left( {\mathit{\boldsymbol{F}}_{\rm{H}}^{ - 1}\mathop {{\mathit{\boldsymbol{E}}_{\rm{H}}}}\limits^\backslash \mathit{\boldsymbol{F}}_{\rm{H}}^{ - {\rm{1}}}} \right) = \mathit{\boldsymbol{\tilde \sigma }}:{\mathit{\boldsymbol{D}}_{\rm{H}}}. \end{array} $ (46)

从式(40) 第二式可知${{\mathit{\boldsymbol{\tilde T}}}_{\rm{H}}}$是EH的热力学共轭量, 式(46) 表明Bishop有效应力${\mathit{\boldsymbol{\tilde \sigma }}}$与孔隙变形率DH成热力学共轭量对.

2 建立耗散率流势函数

把式(39) 代入到式(38) 并利用式(40)~(44) 得

$ \begin{array}{l} \theta \sum\limits_\alpha {{\rho _\alpha }\mathop {{\eta _\alpha }}\limits^\backslash } = \mathit{\boldsymbol{T}}_{\rm{H}}^{\rm{p}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash + {{\tilde s}^{\rm{p}}}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash + \sum\limits_\alpha {P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } + \\ \sum\limits_\alpha {\left( { - {\rm{div}}\;{\mathit{\boldsymbol{q}}_\alpha } + {\rho _\alpha }{r_\alpha }} \right)} + \sum\limits_\beta {{\mathit{\boldsymbol{W}}_\beta } \cdot \left( {{P_\beta }{\rm{grad}}{\varphi _{{\rm{B}}\beta }} - {{\mathit{\boldsymbol{\hat p}}}_\beta }} \right)} . \end{array} $ (47)

利用当Γα=ηα时的式(12), 由式(47) 得

$ \begin{array}{l} \theta \rho \dot \eta = \mathit{\boldsymbol{T}}_{\rm{H}}^{\rm{p}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash + {{\tilde s}^{\rm{p}}}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash + \sum\limits_\alpha {P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } + \\ \sum\limits_\alpha {\left( { - {\rm{div}}\;{\mathit{\boldsymbol{q}}_\alpha } + {\rho _\alpha }{r_\alpha }} \right)} + \sum\limits_\beta {{\mathit{\boldsymbol{W}}_\beta } \cdot \left( {{P_\beta }{\rm{grad}}{\varphi _{{\rm{B}}\beta }} - {{\mathit{\boldsymbol{\hat p}}}_\beta }} \right)} - \\ \theta \sum\limits_\alpha {{\rm{div}}\left( {{\rho _\alpha }{\eta _\alpha }{\mathit{\boldsymbol{u}}_\alpha }} \right)} . \end{array} $ (48)

式中：η为非饱和岩石混合物的熵.令$\boldsymbol{\vartheta} _\theta ^{\rm{P}}$ =θgrad (1/θ), q=$ \sum\limits_\alpha {{q_\alpha }}, \mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}} = {P_\beta }{\rm{grad}}\;{\varphi _{B\beta }} - {{\mathit{\boldsymbol{\hat p}}}_\beta }.$式(48) 变为

$ \begin{array}{l} \rho \dot \eta = \frac{1}{\theta }\left( {\mathit{\boldsymbol{T}}_{\rm{H}}^{\rm{p}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash + {{\tilde s}^{\rm{p}}}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash + \sum\limits_\alpha {P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } + } \right.\\ \;\;\;\;\left. {\sum\limits_\beta {\mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}}} \cdot {\mathit{\boldsymbol{W}}_\beta }} \right) + \frac{1}{\theta }\vartheta _\theta ^{\rm{p}} \cdot \mathit{\boldsymbol{q}} - {\rm{div}}\left( {\frac{\mathit{\boldsymbol{q}}}{\theta } + \sum\limits_\alpha {{\rho _\alpha }{\eta _\alpha }{\mathit{\boldsymbol{u}}_\alpha }} } \right) + \\ \;\;\;\;\sum\limits_\alpha {\frac{{{\rho _\alpha }{r_\alpha }}}{\theta }} . \end{array} $ (49)

式中：α∈{S, L, G};β∈{L, G}.利用式(5) 和(10) 可得[21]

$ \begin{array}{l} \frac{{\partial \left( {\rho \eta } \right)}}{{\partial t}} = \rho \frac{{\partial \eta }}{{\partial t}} + \eta \frac{{\partial \rho }}{{\partial t}} = \\ \;\;\;\;\rho \dot \eta - \rho \mathit{\boldsymbol{v}} \cdot {\rm{grad}}\eta - \eta {\rm{div}}\left( {\rho \mathit{\boldsymbol{v}}} \right) = \rho \dot \eta - {\rm{div}}\left( {\rho \eta \mathit{\boldsymbol{v}}} \right). \end{array} $ (50)

类似于文献[21]的推导, 把式(49) 代入到式(50) 得

$ \begin{array}{l} \frac{{\partial \left( {\rho \eta } \right)}}{{\partial t}} = \underbrace { - {\rm{div}}\left( {\rho \eta \mathit{\boldsymbol{v}} + \sum\limits_\alpha {\frac{{{\mathit{\boldsymbol{q}}_\alpha } + \theta {\rho _\alpha }{\eta _\alpha }{\mathit{\boldsymbol{u}}_\alpha }}}{\theta }} } \right) + \sum\limits_\alpha {\frac{{{\rho _\alpha }{r_\alpha }}}{\theta }} }_{{\rm{entropy}}\;{\rm{flux}}} + \\ \underbrace {\frac{1}{\theta }\left( {\mathit{\boldsymbol{T}}_{\rm{H}}^{\rm{p}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash + {{\tilde s}^{\rm{p}}}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash + \sum\limits_\alpha {P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } } \right) + \sum\limits_\beta {\mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}}} \cdot \frac{{{\mathit{\boldsymbol{W}}_\beta }}}{\theta } + \frac{1}{\theta }\vartheta _\theta ^{\rm{p}} \cdot \mathit{\boldsymbol{q}}}_{{\rm{entropy}}\;{\rm{production}}}. \end{array} $ (51)

从式(51) 可得：非饱和多孔岩石的熵流项为

$ {\eta ^r} = - {\rm{div}}\left( {\rho \eta \mathit{\boldsymbol{v}} + \sum\limits_\alpha {\frac{{{\mathit{\boldsymbol{q}}_\alpha } + \theta {\eta _\alpha }{\rho _\alpha }{\mathit{\boldsymbol{u}}_\alpha }}}{\theta }} } \right) + \sum\limits_\alpha {\frac{{{\rho _\alpha }{r_\alpha }}}{\theta }} . $ (52)

熵产项为

$ \begin{array}{l} {\eta ^i} = \frac{1}{\theta }\left( {\mathit{\boldsymbol{T}}_{\rm{H}}^{\rm{p}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash + {{\tilde s}^{\rm{p}}}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash + \sum\limits_\alpha {P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } } \right) + \\ \;\;\;\;\sum\limits_\beta {\mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}}} \cdot \frac{{{\mathit{\boldsymbol{W}}_\beta }}}{\theta } + \frac{1}{\theta }\vartheta _\theta ^{\rm{p}} \cdot \mathit{\boldsymbol{q}} \ge 0. \end{array} $ (53)

从式(53) 等式的右边可以看出, 第1项(用小括号括出部分)为不可逆变形引起的熵产, 第2项是相对运动导致的力学扩散引起的熵产, 第3项是热扩散引起的熵产.根据式(53) 获得的耗散函数表达式为

$ \begin{array}{l} \mathscr{D} = \theta {\eta ^i} = \left( {\mathit{\boldsymbol{T}}_{\rm{H}}^{\rm{p}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash + {{\tilde s}^{\rm{p}}}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash + \sum\limits_\alpha {P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } } \right) + \\ \;\;\;\;\sum\limits_\beta {P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } + \sum\limits_\beta {\mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}}} \cdot {\mathit{\boldsymbol{W}}_\beta } + \vartheta _\theta ^{\rm{p}} \cdot \mathit{\boldsymbol{q}} = \\ \;\;\;\;\mathit{\boldsymbol{T}}_{\rm{H}}^{\rm{p}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash + {{\tilde s}^{\rm{p}}}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash + \sum\limits_\alpha {P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } + \sum\limits_\beta {\mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}}} \cdot {\mathit{\boldsymbol{W}}_\beta } + \\ \;\;\;\;\vartheta _\theta ^{\rm{p}} \cdot \mathit{\boldsymbol{q}} \ge 0. \end{array} $ (54)

Biot[22]和Housbly等[23]认为存在一个耗散率流势函数$\mathscr{D}^* = \mathscr{D}{^*}(\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash, \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash, \mathop {\vartheta _{_\alpha }^{\rm{p}}}\limits^\backslash, {\mathit{\boldsymbol{W}}_\beta }, \mathit{\boldsymbol{q}})$, 可以把式(54) 中的耗散力表示为

$ \mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}} = \frac{{\partial {\mathscr{D}^ * }}}{{\partial \mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash }},{{\tilde s}^{\rm{p}}} = \frac{{\partial {\mathscr{D}^ * }}}{{\partial \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash }},P_\alpha ^{\rm{p}} = \frac{{\partial {\mathscr{D}^ * }}}{{\partial \mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash }}, $ (55)
$ \mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}} = \frac{{\partial {\mathscr{D}^ * }}}{{\partial {\mathit{\boldsymbol{W}}_\beta }}},\vartheta _\theta ^{\rm{p}} = \frac{{\partial {\mathscr{D}^ * }}}{{\partial \mathit{\boldsymbol{q}}}}. $ (56)

由于$(\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash, \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash, \mathop {\vartheta _{_\alpha }^{\rm{p}}}\limits^\backslash, {\mathit{\boldsymbol{W}}_\beta }, \mathit{\boldsymbol{q}})$是一组与速率有关的热力学流, 故式(55)-(56) 反映了非饱和多孔岩石的不可逆黏性特性.文中的式(40)~(44) 反映的是材料弹性特性以及黏弹性特性的组合关系.式(55)-(56) 揭示的是黏性变形的能耗规律, 若应力一定, 则1) 根据式(40)~(44) 和式(55)-(56) 可获得EH、Sr、${\vartheta _{\rm{S}}}$、${\vartheta _{\rm{L}}}$和${\vartheta _{\rm{G}}}$；2) 再根据CH=2EH+1、CS=CHexp (-2 ${\vartheta _{\rm{S}}}$ /3)、ES=(CS-1)/2、φBS=φBS0/ $\sqrt {{\rm{det}}\left( {{\mathit{\boldsymbol{C}}_{\rm{H}}}} \right)} $获得ES和φBS, 式中的φBS0为固体的初始体积分数；3) 接着根据式(7) 中的φB=1-φBS获得φB, 根据φBL=SrφB和φBG=(1-Sr)φB获得φBL和φBG；4) 再根据ρβ0/ρβ=(φBβ0/φBβ)exp (-${\vartheta _\beta }$)、${\vartheta _\beta }$和φBβ求得ρβ0/ρβ(β∈{L, G}), 式中ρβ0为流体的初始密度, φBβ0为其初始体积分数.鉴于流体的体应变是ρβ0/ρβ的函数, 故根据ρβ0/ρβ可求得液气2个流体的体应变, 这样可全部获得非饱和岩石的3个位移矢量和3个体积分数.式(40)~(44) 与式(55)-(56) 一起反映了非饱和多孔岩石的可逆和不可逆性质, 构成一个基于热力学的非饱和岩石黏弹性本构理论体系.

式(40)~(44) 反映的弹性力学特性还可以采用优化方法表示为

$ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{\xi }} \in \mathit{\boldsymbol{\xi }}_{\rm{S}}^{\rm{e}}} \left[ {{\mathscr{E}_{\rm{S}}} - \frac{1}{{{\rho _{\rm{S}}}}}\left( {\theta {\eta _{\rm{S}}} + {{\mathit{\boldsymbol{\tilde T}}}_{\rm{H}}}:{\mathit{\boldsymbol{E}}_{\rm{H}}} + {\varphi _{{\rm{BS}}}}{P_{\rm{S}}}{\vartheta _{\rm{S}}} - \tilde s{S_{\rm{r}}} + } \right.} \right.\\ \left. {\left. {\mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}}:\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}} + P_{\rm{S}}^{\rm{p}}\vartheta _{\rm{S}}^{\rm{p}} + {{\tilde s}^{\rm{p}}}S_{\rm{r}}^{\rm{p}} + {0_{P_{\rm{L}}^ * - {P_{\rm{L}}}}}{\lambda _{\rm{L}}} + {0_{P_{\rm{G}}^ * - {P_{\rm{G}}}}}{\lambda _{\rm{G}}}} \right)} \right]. \end{array} $ (57)

式中:

$ \mathit{\boldsymbol{\xi }}_{\rm{S}}^{\rm{e}} = \left( {{\eta _{\rm{S}}},{\mathit{\boldsymbol{E}}_{\rm{H}}},{\vartheta _{\rm{S}}},{S_{\rm{r}}},\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}},\vartheta _{\rm{S}}^{\rm{p}},S_{\rm{r}}^{\rm{p}},{\lambda _{\rm{L}}},{\lambda _{\rm{G}}}} \right). $ (58)

和

$ \mathop {\min }\limits_{\mathit{\boldsymbol{\xi }} \in \mathit{\boldsymbol{\xi }}_\beta ^e} \left[ {{\mathit{\mathscr{E}}_\beta } - \frac{1}{{{\rho _\beta }}}\left( {\theta {\eta _\beta } + {\varphi _{{\rm{B}}\beta }}{P_\beta }{\vartheta _\beta } + P_\beta ^{\rm{p}}\vartheta _\beta ^{\rm{p}}} \right)} \right],\beta \in \left\{ {{\rm{L}},{\rm{G}}} \right\}. $ (59)

式中:

$ \mathit{\boldsymbol{\xi }}_\beta ^e = \left( {{\eta _\beta },{\vartheta _\beta },\vartheta _\beta ^{\rm{p}}} \right),\beta = {\rm{L}},{\rm{G}}{\rm{.}} $ (60)

同时, 式(55)-(56) 反映的不可逆力学性质也可以用优化原理表示为

$ \begin{array}{l} \mathop {\max }\limits_{\mathit{\boldsymbol{\xi }} \in {\mathit{\boldsymbol{\xi }}^{\rm{p}}}} \left( {{\mathscr{D}^ * } - \mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash - {{\tilde s}^{\rm{p}}}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash - \sum\limits_\alpha {P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } - } \right.\\ \;\;\;\;\left. {\sum\limits_\beta {\mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}}} \cdot {\mathit{\boldsymbol{W}}_\beta } - \vartheta _\theta ^{\rm{p}} \cdot \mathit{\boldsymbol{q}}} \right). \end{array} $ (61)

式中：

$ {\mathit{\boldsymbol{\xi }}^{\rm{p}}} = \left( {\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash ,\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash ,\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash ,{\mathit{\boldsymbol{W}}_\beta },\mathit{\boldsymbol{q}}} \right). $ (62)

当$\mathscr{D}{^*} =\mathscr{D} {^*}(\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash, \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash, \mathop {\vartheta _{_\alpha }^{\rm{p}}}\limits^\backslash, {\mathit{\boldsymbol{W}}_\beta }, \mathit{\boldsymbol{q}})$是关于($(\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash, \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash, \mathop {\vartheta _{_\alpha }^{\rm{p}}}\limits^\backslash, {\mathit{\boldsymbol{W}}_\beta }, \mathit{\boldsymbol{q}})$的n次齐次函数时, 有

$ \begin{array}{l} n{\mathscr{D}^ * } = \frac{{\partial {\mathscr{D}^ * }}}{{\partial \mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash }}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash + \frac{{\partial {\mathscr{D}^ * }}}{{\partial \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash }}:\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash + \frac{{\partial {\mathscr{D}^ * }}}{{\partial \mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash }}:\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash + \\ \;\;\;\sum\limits_\beta {\frac{{\partial {\mathscr{D}^ * }}}{{\partial {\mathit{\boldsymbol{W}}_\beta }}}} \cdot {\mathit{\boldsymbol{W}}_\beta } + \frac{{\partial {\mathscr{D}^ * }}}{{\partial \mathit{\boldsymbol{q}}}} \cdot \mathit{\boldsymbol{q = }}\\ \;\;\;\mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash + {{\tilde s}^{\rm{p}}}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash + \sum\limits_\alpha {P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } + \sum\limits_\beta {\mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}}} \cdot {\mathit{\boldsymbol{W}}_\beta } + \\ \;\;\;\vartheta _\theta ^{\rm{p}} \cdot \mathit{\boldsymbol{q = }}\mathscr{D} \ge 0. \end{array} $ (63)

当$\mathscr{D}{^*} =\mathscr{D} {^*}(\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash, \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash, \mathop {\vartheta _{_\alpha }^{\rm{p}}}\limits^\backslash, {\mathit{\boldsymbol{W}}_\beta }, \mathit{\boldsymbol{q}})$是关于$(\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash, \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash, \mathop {\vartheta _{_\alpha }^{\rm{p}}}\limits^\backslash, {\mathit{\boldsymbol{W}}_\beta }, \mathit{\boldsymbol{q}})$的一次齐次函数时, 由式(63) 可得

$ \begin{array}{l} {\mathscr{D}^ * } = \mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash + {{\tilde s}^{\rm{p}}}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash + \sum\limits_\alpha {P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } + \sum\limits_\beta {\mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}}} \cdot {\mathit{\boldsymbol{W}}_\beta } + \\ \;\;\;\vartheta _\theta ^{\rm{p}} \cdot \mathit{\boldsymbol{q = }}\mathscr{D} \ge 0. \end{array} $ (64)

式中: $\mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{P}}$, ${{\tilde s}^{\rm{P}}}$, Pαp, $\mathit{\boldsymbol{\hat F}}_\beta ^p$和$\boldsymbol{\vartheta} _\theta ^{\rm{P}}$是EHp, Srp, $\vartheta _\alpha ^{\rm{P}}$, Wβ和q关于耗散率流势函数的热力学共轭量, 故相应地耗散率力势函数可表示为

$ \begin{array}{l} {f^ * } = \left( {\mathit{\boldsymbol{T}}_{\rm{H}}^{\rm{p}},{{\tilde s}^{\rm{p}}},P_\alpha ^{\rm{p}},\mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}},\vartheta _\theta ^{\rm{p}}} \right) = {\mathscr{D}^ * } - \frac{{\partial {\mathscr{D}^ * }}}{{\partial \mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash }}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash - \\ \;\;\;\frac{{\partial {\mathscr{D}^ * }}}{{\partial \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash }}:\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash - \sum\limits_\alpha {\frac{{\partial {\mathscr{D}^ * }}}{{\partial \mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash }}:\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } - \sum\limits_\beta {\frac{{\partial {\mathscr{D}^ * }}}{{\partial {\mathit{\boldsymbol{W}}_\beta }}}} \cdot {\mathit{\boldsymbol{W}}_\beta } - \\ \;\;\;\frac{{\partial {\mathscr{D}^ * }}}{{\partial \mathit{\boldsymbol{q}}}} \cdot \mathit{\boldsymbol{q}} = {\mathscr{D}^ * } - \mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash - {{\tilde s}^{\rm{p}}}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash - \sum\limits_\alpha {P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } - \\ \;\;\;\sum\limits_\beta {\mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}}} \cdot {\mathit{\boldsymbol{W}}_\beta } - \vartheta _\theta ^{\rm{p}} \cdot \mathit{\boldsymbol{q = }}0. \end{array} $ (65)

故当$\mathscr{D} {^*}(\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash, \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash, \mathop {\vartheta _{_\alpha }^{\rm{p}}}\limits^\backslash, {\mathit{\boldsymbol{W}}_\beta }, \mathit{\boldsymbol{q}})$是关于$(\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash, \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash, \mathop {\vartheta _{_\alpha }^{\rm{p}}}\limits^\backslash, {\mathit{\boldsymbol{W}}_\beta }, \mathit{\boldsymbol{q}})$一次齐次函数时, 由式(61)-(62) 结合式(65) 可得

$ \begin{array}{l} \mathop {\max }\limits_{\mathit{\boldsymbol{\xi }} \in {\mathit{\boldsymbol{\xi }}^{\rm{p}}}} \left( {{f^ * } - \mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash :\mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}} - \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash {{\tilde s}^{\rm{p}}} - \sum\limits_\alpha {\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash P_\alpha ^{\rm{p}}} - } \right.\\ \;\;\;\;\left. {\sum\limits_\beta {{\mathit{\boldsymbol{W}}_\beta } \cdot \mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}} - \mathit{\boldsymbol{q}} \cdot \vartheta _\theta ^{\rm{p}}} } \right)\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;{f^ * } = \left( {\mathit{\boldsymbol{T}}_{\rm{H}}^{\rm{p}},{{\tilde s}^{\rm{p}}},P_\alpha ^{\rm{p}},\mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}},\vartheta _\theta ^{\rm{p}}} \right) = 0. \end{array} $ (66)

式(66) 类似于Hill最大耗散功原理[24].由式(66) 可得

$ \mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash = \mathop {{\lambda ^{\rm{p}}}}\limits^\backslash \frac{{\partial {f^ * }}}{{\partial \mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}}}},\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash = \mathop {{\lambda ^{\rm{p}}}}\limits^\backslash \frac{{\partial {f^ * }}}{{\partial {{\tilde s}^{\rm{p}}}}},\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash = \mathop {{\lambda ^{\rm{p}}}}\limits^\backslash \frac{{\partial {f^ * }}}{{\partial P_\alpha ^{\rm{p}}}}. $ (67)
$ {\mathit{\boldsymbol{W}}_\beta } = \mathop {{\lambda ^{\rm{p}}}}\limits^\backslash \frac{{\partial {f^ * }}}{{\partial \mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}}}},\mathit{\boldsymbol{q = }}\mathop {{\lambda ^{\rm{p}}}}\limits^\backslash \frac{{\partial {f^ * }}}{{\partial \vartheta _\theta ^{\rm{p}}}}. $ (68)

式中:λp为塑性因子.

当变量($\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash, \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash, \mathop {\vartheta _{_{\rm{S}}}^{\rm{p}}}\limits^\backslash $)与($\mathop {\vartheta _{_\beta }^{\rm{p}}}\limits^\backslash, {\mathit{\boldsymbol{W}}_\beta }, \mathit{\boldsymbol{q}}$)相互独立时(β∈{L, G}), 由式(54) 可得

$ {\mathscr{D}_{\rm{p}}} = \mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash + {{\tilde s}^{\rm{p}}}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash + P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash \ge 0, $ (69)
$ {\mathscr{D}_{\rm{m}}} = \sum\limits_\alpha {P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash } + \sum\limits_\beta {{\mathit{\boldsymbol{W}}_\beta } \cdot \mathit{\boldsymbol{\hat F}}_\beta ^{\rm{p}} + \mathit{\boldsymbol{q}} \cdot \vartheta _\theta ^{\rm{p}}} \ge 0. $ (70)

式(54) 中的$\mathscr{D} = {\mathscr{D}_{\rm{P}}} + {\mathscr{D}_{\rm{m}}}$.由此Biot建议的耗散率流势函数可表示为$\mathscr{D}{^*} = \mathscr{D}{_{\rm{p}}}^*(\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash, \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash, \mathop {\vartheta _{_{\rm{S}}}^{\rm{p}}}\limits^\backslash ) + \mathscr{D}{_{\rm{m}}}(\mathop {\vartheta _{_\beta }^{\rm{p}}}\limits^\backslash, {\mathit{\boldsymbol{W}}_\beta }, \mathit{\boldsymbol{q}})$.若$\mathscr{D}{_{\rm{p}}}^*(\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash, \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash, \mathop {\vartheta _{_{\rm{S}}}^{\rm{p}}}\limits^\backslash )$是关于($\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash, \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash, \mathop {\vartheta _{_{\rm{S}}}^{\rm{p}}}\limits^\backslash $)的一次齐次欧拉函数, 则与式(64)~(68) 的推导类似, 令

$ \begin{array}{l} f_{\rm{p}}^ * = \left( {\mathit{\boldsymbol{T}}_{\rm{H}}^{\rm{p}},{{\tilde s}^{\rm{p}}},P_{\rm{S}}^{\rm{p}}} \right) = \mathscr{D}_{\rm{p}}^ * - \frac{{\partial \mathscr{D}_{\rm{p}}^ * }}{{\partial \mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash }}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash - \frac{{\partial {\mathscr{D}^ * }}}{{\partial \mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash }}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash - \\ \;\;\;\;\;\frac{{\partial \mathscr{D}_{\rm{p}}^ * }}{{\partial \mathop {\vartheta _{\rm{S}}^{\rm{p}}}\limits^\backslash }}\mathop {\vartheta _{\rm{S}}^{\rm{p}}}\limits^\backslash = \mathscr{D}_{\rm{p}}^ * - \mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}}:\mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash - {{\tilde s}^{\rm{p}}}\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash - P_\alpha ^{\rm{p}}\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash = 0. \end{array} $ (71)

由此可得

$ \mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash = \mathop {{\lambda ^{\rm{p}}}}\limits^\backslash \frac{{\partial f_{\rm{p}}^ * }}{{\partial \mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}}}},\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash = \mathop {{\lambda ^{\rm{p}}}}\limits^\backslash \frac{{\partial f_{\rm{p}}^ * }}{{\partial {{\tilde s}^{\rm{p}}}}},\mathop {\vartheta _\alpha ^{\rm{p}}}\limits^\backslash = \mathop {{\lambda ^{\rm{p}}}}\limits^\backslash \frac{{\partial f_{\rm{p}}^ * }}{{\partial P_{\rm{S}}^{\rm{p}}}}. $ (72)

式(71)-(72) 揭示了塑性位势理论的热力学依据.当fp*进一步与$\vartheta _{\rm{S}}^{\rm{p}}$也无关时, 则式(72) 变为

$ \mathop {\mathit{\boldsymbol{E}}_{\rm{H}}^{\rm{p}}}\limits^\backslash = \mathop {{\lambda ^{\rm{p}}}}\limits^\backslash \frac{{\partial f_{\rm{p}}^ * }}{{\partial \mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}}}},\mathop {S_{\rm{r}}^{\rm{p}}}\limits^\backslash = \mathop {{\lambda ^{\rm{p}}}}\limits^\backslash \frac{{\partial f_{\rm{p}}^ * }}{{\partial {{\tilde s}^{\rm{p}}}}}. $ (73)

式(73) 与经典塑性力学位势理论相同.

现在假定混合物均匀化响应原理成立, 根据式(34) 和(38) 之间的文字分析, 此时孔隙和岩石材料的变形相互解耦, 则式(39) 中多孔固相的自由能可分为孔隙自由能$\mathscr{E}_{\text{SH}}$和材料自由能$\mathscr{E}_{\text{SM}}$ 2个部分.多孔固相自由能$\mathscr{E}_{\text{S}}$ (ηS, EH, Sr, ${\vartheta _{\rm{S}}}$, EHp, Srp, $\vartheta _{\rm{S}}^{\rm{p}}$, λL, λG)可表示为

$ \begin{array}{l} {\mathscr{E}_{\rm{S}}} = {\mathscr{E}_{{\rm{SH}}}}\left( {{\eta _{{\rm{SH}}}},{E_{\rm{H}}},{S_{\rm{r}}},E_{\rm{H}}^{\rm{p}},S_{\rm{r}}^{\rm{p}},{\lambda _{\rm{L}}},{\lambda _{\rm{G}}}} \right) + \\ \;\;\;\;\;\;\;\;{\mathscr{E}_{{\rm{SM}}}}\left( {{\eta _{{\rm{SM}}}},{\vartheta _{\rm{S}}},\vartheta _{\rm{S}}^{\rm{p}}} \right). \end{array} $ (74)

式中:ηS=ηSH+ηSM.把式(74) 代入到式(40)~(44) 得

$ \theta = \frac{{\partial {\mathscr{E}_{{\rm{SH}}}}}}{{\partial {\eta _{{\rm{SH}}}}}} = \frac{{\partial {\mathscr{E}_{{\rm{SM}}}}}}{{\partial {\eta _{{\rm{SM}}}}}},{{\mathit{\boldsymbol{\tilde T}}}_{\rm{H}}} = {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{{\rm{SH}}}}}}{{\partial {\mathit{E}_{\rm{H}}}}},\tilde s = - {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{{\rm{SH}}}}}}{{\partial {S_{\rm{r}}}}}, $ (75)
$ \mathit{\boldsymbol{\tilde T}}_{\rm{H}}^{\rm{p}} = - {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{{\rm{SH}}}}}}{{\partial \mathit{E}_{\rm{H}}^{\rm{p}}}},\;\;\;\;\;{{\tilde s}^{\rm{p}}} = - {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{{\rm{SH}}}}}}{{\partial S_{\rm{r}}^{\rm{p}}}}, $ (76)
$ {0_{P_{\rm{L}}^ * - {P_{\rm{L}}}}} = {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{{\rm{SH}}}}}}{{\partial {\lambda _{\rm{L}}}}},{0_{P_{\rm{G}}^ * - {P_{\rm{G}}}}} = {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{{\rm{SH}}}}}}{{\partial {\lambda _{\rm{G}}}}}, $ (77)
$ {P_{\rm{S}}} = {\rho ^{\rm{S}}}\frac{{\partial {\mathscr{E}_{{\rm{SM}}}}}}{{\partial {\vartheta _{\rm{S}}}}},P_{\rm{S}}^{\rm{p}} = - {\rho _{\rm{S}}}\frac{{\partial {\mathscr{E}_{{\rm{SM}}}}}}{{\partial \vartheta _{\rm{S}}^{\rm{p}}}}. $ (78)

式(74)~(78) 表明, 当岩石的孔隙变形与材料变形相互独立时, 固相自由能可表示为改变岩石裂隙变化的自由能和产生岩石材料变形的自由能之和, 式(74)~(76) 表明裂隙变化自由能取决于孔隙Green应变、孔隙Green内应变、饱和度和内饱和度, 式(74) 和(78) 表明岩石材料自由能取决于岩石的材料体应变和材料内体应变.

Borja主张式(31) 中等式右边的前2项满足

$ \mathit{\boldsymbol{\tilde \sigma }}:{\rm{grad}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} + \bar P{\varphi _{{\rm{BS}}}}\frac{{\mathop {{\rho ^{\rm{S}}}}\limits^\backslash }}{{{\rho ^{\rm{S}}}}} = \left[ {\mathit{\boldsymbol{\sigma + }}\left( {1 - \frac{{{K_{\rm{b}}}}}{{{K_{\rm{S}}}}}} \right)\bar P{\bf{1}}} \right]:{\rm{grad}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}}. $ (79)

式中：Kb为裂隙岩石的压缩模量, KS为岩石材料的压缩模量, 根据上式、${\vartheta _{\rm{S}}}$的定义和$\mathit{\boldsymbol{\tilde \sigma = \sigma + }}\bar P{\bf{1}}$可得

$ {\varphi _{{\rm{BS}}}}\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash = - \frac{{{K_{\rm{b}}}}}{{{K_{\rm{S}}}}}{\rm{grad}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}}:{\bf{1}} = - \frac{{{K_{\rm{b}}}}}{{{K_{\rm{S}}}}}{\rm{div}}\left( {{\mathit{\boldsymbol{v}}_{\rm{S}}}} \right). $ (80)

可惜的是, 式(80) 在绝大多数条件下不成立.现按最简单的情况考虑：混合物均匀化响应原理成立且非饱和裂隙岩石发生小应变线弹性变形且不考虑饱和度对裂隙岩石变形的影响.令

$ \begin{array}{l} {\mathscr{E}_{\rm{S}}} = \frac{{{K_{\rm{H}}}}}{{2{\rho _{\rm{S}}}}}{\left( {{\mathit{\boldsymbol{E}}_{\rm{H}}}:{\bf{1}}} \right)^2} + \frac{{P_{\rm{L}}^ * - {P_{\rm{L}}}}}{{{\rho _{\rm{S}}}}}{\lambda _{\rm{L}}} + \frac{{P_{\rm{G}}^ * - {P_{\rm{G}}}}}{{{\rho _{\rm{S}}}}}{\lambda _{\rm{G}}} + \\ \;\;\;\;\;\;\;\;\frac{{{K_{\rm{S}}}}}{{2{\rho ^{\rm{S}}}}}\vartheta _{\rm{S}}^2. \end{array} $ (81)

式中：KH为孔隙体积变形, 根据小应变条件(此时ρS和ρS可视为常量, ${{\tilde \sigma }_{\rm{m}}} \simeq {{\mathit{\boldsymbol{\tilde T}}}_{\rm{H}}}:{\bf{1}}/3, {\rm{div}}\;\;\;\mathit{\boldsymbol{v}}_{\text{s}} + \mathop {{\vartheta _{\rm{S}}}}\limits^\backslash \simeq {\mathit{\boldsymbol{E}}_{\rm{H}}}{\bf{1}})$、式(74)-(75)、(78) 和上述条件可容易地获得

$ {\rm{div}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} + \mathop {{\vartheta _{\rm{S}}}}\limits^\backslash = \frac{1}{{{K_{\rm{H}}}}}\mathop {{{\mathit{\boldsymbol{\tilde \sigma }}}_{\rm{m}}}}\limits^\backslash ,\;\;\;\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash = \frac{1}{{{K_{\rm{S}}}}}\mathop {{P_{\rm{S}}}}\limits^\backslash . $ (82)

根据式(77) 可知式(35) 进而式(27) 成立.根据Kb是在平均孔压等于零的条件下测定的, 有

$ \frac{1}{{{K_{\rm{H}}}}} = \frac{1}{{{K_{\rm{b}}}}} - \frac{1}{{{\varphi _{{\rm{BS}}}}{K_{\rm{S}}}}}. $ (83)

注意到小应变条件下φBS和φBF可视为常数, 把式(83) 代入到式(82) 并利用${P_{\rm{T}}} = \bar P - {{\tilde \sigma }_{\text{m}}}$得

$ \begin{array}{l} {\rm{div}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} = \left( {\frac{1}{{{K_{\rm{b}}}}} - \frac{1}{{{\varphi _{{\rm{BS}}}}{K_{\rm{S}}}}}} \right)\mathop {{{\tilde \sigma }_{\rm{m}}}}\limits^\backslash - \mathop {{\vartheta _{\rm{S}}}}\limits^\backslash = \\ \;\;\;\left( {\frac{1}{{{K_{\rm{b}}}}} - \frac{1}{{{\varphi _{{\rm{BS}}}}{K_{\rm{S}}}}}} \right)\left( {\mathop {\bar P}\limits^\backslash - \mathop {{P_{\rm{T}}}}\limits^\backslash } \right) - \frac{1}{{{\varphi _{{\rm{BS}}}}{K_{\rm{S}}}}}\left( {\mathop {{P_{\rm{T}}}}\limits^\backslash - {\varphi _{\rm{B}}}\mathop {\bar P}\limits^\backslash } \right) = \\ \;\;\; - \frac{1}{{{K_{\rm{b}}}}}\left[ {\mathop {{P_{\rm{T}}}}\limits^\backslash - \left( {1 - \frac{{{K_{\rm{b}}}}}{{{K_{\rm{S}}}}}} \right)\mathop {\bar P}\limits^\backslash } \right]. \end{array} $ (84)

式(84) 等同于非饱和多孔岩石的Skempton有效应力公式(陈正汉非饱和多孔介质有效应力公式[14]或Borja非饱和多孔介质有效应力公式[15]), 但即使在这样的条件下下列不等式依然成立：

$ \begin{array}{l} {\varphi _{{\rm{BS}}}}\mathop {{\vartheta _{\rm{S}}}}\limits^\backslash = \frac{{\left( {\mathop {{P_{\rm{T}}}}\limits^\backslash - {\varphi _{\rm{B}}}\mathop {\bar P}\limits^\backslash } \right)}}{{{K_{\rm{S}}}}} \ne \frac{{{K_{\rm{b}}}}}{{{K_{\rm{S}}}}}\frac{1}{{{K_{\rm{b}}}}}\left[ {\mathop {{P_{\rm{T}}}}\limits^\backslash - } \right.\\ \left. {\left( {1 - \frac{{{K_{\rm{b}}}}}{{{K_{\rm{S}}}}}} \right)\mathop {\bar P}\limits^\backslash } \right] = - \frac{{{K_{\rm{b}}}}}{{{K_{\rm{S}}}}}{\rm{div}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}} = - \frac{{{K_{\rm{b}}}}}{{{K_{\rm{S}}}}}{\rm{grad}}\;{\mathit{\boldsymbol{v}}_{\rm{S}}}:{\bf{1}}. \end{array} $ (85)

式(85) 不等式成立, 故式(80) 进而式(79) 不成立.这一结论表明即使类似Skempton公式的非饱和岩土有效应力公式[15]成立, 但非饱和多孔介质的能量平衡方程式(30) 依然不能用类似Skempton公式的非饱和岩土有效应力公式来表达.

3 结论

本文应用混合物理论和不可逆过程热力学, 建立了一个能够同时考虑多孔岩石和岩石材料不可逆变形的非饱和多孔岩石本构理论框架, 获得以下的研究成果：

(1) 根据传统混合物理论, 通过质量平衡方程和组分体积分数, 获得采用Bishop有效应力和组分材料球应力表示的非饱和多孔岩石的能量平衡方程, 指出决定非饱和裂隙岩石本构特性的状态应变量是孔隙应变、饱和度、岩石的材料体应变和液气2相的材料体应变, 相应的状态应力量是Bishop有效应力、有效吸力、岩石的材料球应力和液气2相的材料球应力(孔压), 由此建立了非饱和多孔岩石的自由能势函数.

(2) 确定了非饱和多孔岩石的熵流和熵产, 根据熵产表达式的热力学流和热力学力, 建立了非饱和多孔岩石的耗散率势函数.与自由能势函数一起, 构成非饱和岩石的超黏弹性本构理论框架.

(3) 利用耗散率势函数和一次欧拉函数的性质, 证明本文非饱和多孔岩石超黏性理论包含了基于塑性势的非饱和多孔岩石塑性理论, 为今后建立能够同时考虑多孔岩体和岩石材料弹性和塑性变形的具体本构模型打下坚实的理论基础.

参考文献
[1] LI X S. Thermodynamics-based constitutive frame-work for unsaturated soils[J]. Géotechnique, 2007, 57(5): 411–422. DOI:10.1680/geot.2007.57.5.411
[2] 缪林昌. 非饱和土的本构模型研究[J]. 岩土力学, 2007, 28(5): 855–860.
MIAO Lin-chang. Research of constitutive model of unsaturated soils[J]. Rock and Soil Mechanics, 2007, 28(5): 855–860.
[3] 孙德安, 高游. 不同制样方法非饱和土的持水特性研究[J]. 岩土工程学报, 2015, 37(1): 91–97.
SUN De-an, GAO You. Water retention behaviour of soils with different preparations[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(1): 91–97. DOI:10.11779/CJGE201501010
[4] MA T T, WEI C F, WEI H Z, et al. Hydraulic and mechanical behavior of unsaturated silt: experimental and theoretical characterization[J]. International Journal of Geomechanics, 2015: D4015007-1–D4015007-13.
[5] 胡亚元. 考虑吸附水的非饱和土耗散本构关系研究[J]. 岩土力学, 2016, 36(S1): 14–18.
HU Ya-yuan. Study on the dissipative constitutive relation of unsaturated soil considering adsorbed water[J]. Rock and Soil Mechanics, 2016, 36(S1): 14–18.
[6] 韩同春, 豆红强, 马世国, 等. 考虑雨水重分布对均质无限长边坡稳定性的研究[J]. 浙江大学学报:工学版, 2013, 47(10): 1824–1829.
HAN Tong-chun, DOU Hong-qiang, MA Shi-guo, et al. Rainwater redistribution on stability of homogen-ous infinite slop[J]. Journal of Zhejiang University: Engineering Science, 2013, 47(10): 1824–1829.
[7] 周创兵, 陈益峰, 姜清辉, 等. 复杂岩体多场广义耦合分析导论[M]. 北京: 中国水利水电出版社, 2008: 152-187-328-363.
[8] 赵阳升. 多孔介质多场耦合作用及其工程响应[M]. 北京: 科学出版社, 2010: 171-310.
[9] SKEMPTON A W. Effective stress in soils, concrete and rocks[M]. London, UK: Butterwoths, 1961: 4-16.
[10] SUKLJE L. Rheological aspects of soil mechanics[M]. Interscience, New York, 1969: 123.
[11] DE BUHAM P, DORMIEUX L. On the validity of the effective stress concept for assessing the strength of saturated porous materials: A homogenization approach[J]. Journal of the Mechanics and Physics of Solids, 1996, 44(10): 1649–1667. DOI:10.1016/0022-5096(96)00046-4
[12] WALSH J, BROWN S, DURHAM W. Effective media theory with spatial correlation for flow in fracture[J]. Journal. Geophysical Resiews, 1997, 102(B10): 22587–22594. DOI:10.1029/97JB01895
[13] CARIOU S, DUAN Z, DAVY C, et al. Poro-mechanics of partially saturated COx argillite[J]. Applied Clay Science, 2012, 56: 36–47. DOI:10.1016/j.clay.2011.11.021
[14] 陈正汉. 非饱和土与特殊土力学的基本理论研究[J]. 岩土工程学报, 2014, 36(2): 201–272.
CHEN Zheng-han. On basic theories of unsaturated soils and special soils[J]. Chinese Journal of Geotechni-cal Engineering, 2014, 36(2): 201–272. DOI:10.11779/CJGE201402001
[15] BORJA R I. On the mechanical energy and effective stress in saturated and unsaturated porous continua[J]. International Journal of Solids and Structures, 2006, 43: 1764–1786. DOI:10.1016/j.ijsolstr.2005.04.045
[16] BISHOP A W. The principle of effective stress[J]. Teknisk Ukeblad, 1959, 106(39): 113–143.
[17] 陈卫忠, 谭贤君, 伍国军, 等. 非饱和岩石温度-渗流-应力耦合模型研究[J]. 岩石力学与工程学报, 2007, 26(12): 2395–2403.
CHEN Wei-zhong, TAN Xian-jun, WU Guo-jun, et al. Study on thermo-hydro-mechanical coupling model for unsaturated rock[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(12): 2395–2403. DOI:10.3321/j.issn:1000-6915.2007.12.003
[18] 赵成刚, 刘艳. 连续孔隙介质土力学及其在非饱和土本构关系中的应用[J]. 岩土工程学报, 2009, 31(9): 1324–1335.
ZHAO Cheng-Gang, LIU Yan. Continuum porous medium soil mechanics and its application in constitutive relationship of unsaturated soils[J]. Chinese Journal of Geotechnical Engineering, 2009, 31(9): 1324–1335.
[19] BOWEN R M著, 徐慧己, 张志新, 李如庆等译. 混合物理论[M]. 南京: 江苏科学技术出版社, 1983: 1-48.
[20] GEERTSMA J. The effect of fluid pressure decline on volumetric changes of porous rocks[J]. Petroleum Transactions, 1957, 210: 331–339.
[21] 李如生. 非平衡态热力学和耗散结构[M]. 北京: 清华大学出版社, 1986: 76-106.
[22] BO IT, M .A. Theory of stress-strain relations in anisotropic viscoelasticity and relaxation phenomena[J]. Journal of Applied Physics, 1954, 25: 1385–1391. DOI:10.1063/1.1721573
[23] HOULSBY G T, PUZRIN A M. Rate-dependent plasticity models derived from potential functions[J]. Journal of Rheology, 2002, 46(1): 113–126. DOI:10.1122/1.1427911
[24] HILL R. The mathematical theory of plasticity[M]. [S.l.]: Oxford University Press, 1950, 50-52.