0.   引言

在核废料处置库、垃圾填埋场、石油和天然气储油层以及海洋工程等涉及环境和岩土相互作用的工程建设领域，经常会遇到盐溶液饱和岩土材料^[1-13]。与普通饱和岩土相比，盐溶液饱和岩土中的可溶盐分在溶液中分解出阴阳离子。这些离子具有化学活性，能够降低矿物表面双电层的厚度，改变岩土颗粒间的距离^[7]，进而引起岩土材料渗透、变形和强度性质发生变化，诱发其工程性能劣化甚至失效。因此研究盐溶液饱和岩土材料的化学活性和工程力学特性的耦合作用具有重要的理论和工程实践意义。

目前国内外学者已经较系统地开展了盐溶液饱和岩土材料工程特性的试验研究^[4-6]，并已开始着手盐溶液饱和岩土材料的本构建模研究工作^[8-13]。Guimarães等^[8]基于唯象理论，采用盐溶液中离子摩尔浓度和渗透压作为化学状态变量来建立非饱和黏土的本构模型。徐永福^[3]用表面分形学来完善双电层理论，利用渗透压修正有效应力，以此计算NaCl溶液饱和膨润土的膨胀变形。Dominijanni等^[9]基于固相骨架电荷与孔隙溶液离子之间的电荷平衡条件，引入渗透压作为化学状态变量来建立理论模型。Bennethum等^[10]在Hassanizadeh等^[11]所创立的复合混合物理论基础之上，定义了与经典热力学相容的化学势，通过应变在各组分自由能之间的耦合作用来反映多孔介质的吸附作用，建立了膨胀多孔骨架的本构理论框架。Wei^[2]和Ma等^[12]从混合物理论出发，推导了考虑表面势的粒间应力公式，进而利用渗透压来修正有效应力，建立了考虑化学活性影响的白垩石吸水本构模型。上述成果有力地推动了饱和与非饱和岩土材料在力学-渗流-化学多场耦合研究领域的发展。

在当前考虑化学活性的研究中，大多数学者选用渗透压作为化学状态变量来反映饱和岩土材料中的水-力-化学之间的耦合作用。然而，渗透压的物理含义是半透膜达到热力学平衡时两侧存在的压力差^[3]。但在饱和岩土固结过程中，孔隙流体存在渗流速度，因而未达到热力学平衡状态，与渗透压的物理含义略有不符。除此之外，也有学者如Loret等^[13]基于热力学理论，选用组分的化学势与摩尔分数作为化学状态变量建立了饱和黏土的化学-力学耦合本构模型。但在该理论中，摩尔分数定义在初始构形的摩尔浓度之上，而在试验中能够直接测量的是现时构形的摩尔浓度，这给试验测量带来一定的困难。因此，需要基于现时构形，开展考虑力学-渗流-化学多场耦合的理论建模研究工作。

另外，作为多孔介质，饱和岩土中的孔隙率反映了颗粒组成结构，在水-力-化学耦合机制中起到至关重要的作用^[14]。然而，对于岩土材料中的裂隙岩体，岩石材料本身的变形不可忽略^[15-16]，故对于一般的岩土材料，固相应变包含了固相材料变形引起的基质应变和孔隙变化引起的骨架应变^[15-16]，它们具有不同的力学性质和本构关系。因此需要把基质应变和骨架应变从固相应变中单独区分开来，以便揭示孔隙结构在多孔介质水-力-化学作用过程中的关键作用。此外，当前盐溶液饱和岩土材料的本构研究多侧重于固相组分的本构模型，而忽略了流相组分的本构模型和溶质渗流-扩散控制方程。事实上，饱和岩土的本构关系不仅应包含固相的本构关系，还应包含流相的本构关系以及溶质的渗流-扩散控制方程^[17-19]。由于变形过程中流固两相的体积分数之和必须等于1，同时建立固相和流相的本构模型可以很容易发现它们均包含骨架体积应变，因此更有利于直观揭示流固两相的几何耦合机制^{[15-16, 18-19]}。

在以往的饱和岩土材料本构模型研究历史中，Collins等^[20]、Borja^[21]、陈正汉^[22-23]、赵成刚等^[24]、刘艳等^[25]和胡亚元^[15-16]在混合物理论与非平衡态热力学理论框架下，从热力学第一定律所蕴含的自由能和热力学第二定律所蕴含的耗散势中出发，成功地建立了考虑力学-渗流耦合的饱和岩土弹塑性本构理论。其中自由能决定材料的弹性性质，耗散势则决定材料的塑性性质。本文将在上述研究基础上，把热力学溶液理论^[26]引入混合物理论与非平衡态热力学理论框架，建立考虑力学-渗流-化学多场耦合过程的盐溶液饱和多孔介质的本构理论。

1.   体积分数和密度

盐溶液饱和多孔介质包括固相与孔隙流相两相，为简便分析，把固相视为单一组分构成，流相视为溶质与溶剂两种组分构成。由于可能存在溶解（或析出）、吸附和化学反应等物理化学作用，为不失一般性，假定固相与流相之间以及流相中各组分之间存在质量交换。令α={S, F}为相的指征变量，固相用S表示，流相用F表示；${\varphi _\alpha }$表示$\alpha$相的现时体积分数，${\rho _\alpha }$为$\alpha$相的表观密度^[27]（有些文献也称为组分密度），${\rho _{{\text{R}}\alpha }}$为$\alpha$相的真实密度（或称材料密度）且满足${\rho _\alpha } = {\varphi _\alpha }{\rho _{{\text{R}}\alpha }}$。则饱和多孔介质的总密度为$\rho = {\rho _{\text{S}}} + {\rho _{\text{F}}}$。根据体积分数的定义有${\varphi _{\text{S}}} + {\varphi _{\text{F}}} = 1$。

孔隙流相中包含溶质与溶剂两种组分，根据经典热力学溶液理论^[26]可知，溶质与溶剂相互融合，很难区分出其各自所占体积，因此可认为溶质与溶剂具有与流相相同的体积。令β={L, c}为流相中组分的指征变量，溶剂用L表示，溶质用c表示，$\rho _{\text{F}}^\beta$为流相中β组分的表观密度^[27]，根据表观密度的定义有

2.   平衡方程

2.1   质量守恒

根据混合物理论^[17]可知，盐溶液饱和多孔介质的固相和流相连续，并共享混合物的空间位置。令$\alpha$相的初始位置为${{\boldsymbol{X}}_\alpha }$，$t$时刻的空间位置为${\boldsymbol{x}}$，则$\alpha$相的运动方程可表示为${\boldsymbol{x}} = {{\boldsymbol{x}}_\alpha }({{\boldsymbol{X}}_\alpha }, t)$，令${{\boldsymbol{v}}_\alpha }$为$\alpha$相的速度，则饱和多孔介质混合物的相平均速度为${\boldsymbol{v}} = ({\rho _{\text{S}}}{{\boldsymbol{v}}_{\text{S}}} + {\rho _{\text{F}}}{{\boldsymbol{v}}_{\text{F}}})/\rho$；$\alpha$相相对于多孔介质混合物的相对速度为${{\boldsymbol{u}}_\alpha } = {{\boldsymbol{v}}_\alpha } - {\boldsymbol{v}}$。选择固相作为盐溶液饱和多孔介质的参考构形，令${{\boldsymbol{W}}_{\text{F}}} = {{\boldsymbol{\nu }}_{\text{F}}} - {{\boldsymbol{\nu }}_{\text{S}}}$为流相相对于固相的扩散速度。考虑到可能存在的溶解（或析出）和化学反应等物理化学作用，固流两相之间以及溶剂与溶质之间可能存在质量交换，令${\hat c_{\text{S}}}$与${\hat c_{\text{F}}}$分别为固相和流相之间的质量交换率，$\hat c_{\text{F}}^{\text{L}}$与$\hat c_{\text{F}}^{\text{c}}$分别为流相中溶剂、溶质与固相组分之间的质量交换率。定义质量交换率以质量供给为正。根据以往研究，在考虑固流两相质量交换情况下可得固流两相的质量守恒方程为^{[15, 17]}

式（2）满足${\hat c_{\text{S}}} + {\hat c_{\text{F}}} = 0$^[17]，同理，可得流相中各组分的质量守恒方程为^[10]

式中：$\hat c_{\text{F}}^{\text{L}} + \hat c_{\text{F}}^{\text{c}} = {\hat c_{\text{F}}}$^[10]；$c_{\text{F}}^\beta = \rho _{\text{F}}^\beta /{\rho _{\text{F}}}$为流相中β组分关于现时构形的质量分数；${\boldsymbol{v}}_{\text{F}}^\beta$为流相中β组分的速度；${\boldsymbol{u}}_{\text{F}}^\beta = {\boldsymbol{v}}_{\text{F}}^\beta - {{\boldsymbol{v}}_{\text{F}}}$为流相中β组分相对于流相的相对速度。

令${{\boldsymbol{\sigma }}_\alpha }(\alpha \in \{ {\text{S, F}}\} )$为$\alpha$相的表观应力（有些文献也称为组分应力），${\boldsymbol{\sigma }}$为总应力，为了更好地应用于岩土工程，根据岩土力学规定，应力张量和应变张量均以受压为正。根据混合物理论^[17]与复杂介质的平均应力定理^[28]可得

令总压力为$P_{\mathrm{T}}=\boldsymbol{\sigma}: \boldsymbol{I} / 3$，$P_{\mathrm{S}}=\boldsymbol{\sigma}_{\mathrm{S}}: \boldsymbol{I} /\left(3 \varphi_{\mathrm{S}}\right)$为固相基质压力；$P_{\mathrm{F}}=\boldsymbol{\sigma}_{\mathrm{F}}: \boldsymbol{I} /\left(3 \varphi_{\mathrm{F}}\right)$为流相基质压力或称为孔压。则根据式（4）有${P_{\text{T}}} = {\varphi _{\text{S}}}{P_{\text{S}}} + {\varphi _{\text{F}}}{P_{\text{F}}}$。

2.2   能量平衡

令$\boldsymbol{q}_\alpha$，${r_\alpha }$以及${\hat \varepsilon _\alpha }$分别为$\alpha$相的热流向量、外热供给量以及能量供给量，${\xi _\alpha }$为内能密度；${\boldsymbol{b}}_{\text{F}}^\beta$为外体力密度，则固相与流相的能量平衡方程为^[15-17]

令$\boldsymbol{D}_{\mathrm{S}}=\left[\nabla \boldsymbol{v}_{\mathrm{S}}+\left(\nabla \boldsymbol{v}_{\mathrm{S}}\right)^{\mathrm{T}}\right] / 2$为固相变形率；根据混合物理论^[17]可知，式（5）中的${\hat \varepsilon _\alpha }$满足方程：

孔隙率在多孔介质中反映微观结构的变化，在多孔介质变形中占据着关键地位^[14]。在土力学中，把与孔隙率变化存在一一对应关系的体积应变称为骨架应变^[15-16]。因此，有必要将孔隙率变化相关的骨架应变从固相应变中分离开来。然而下面的分析将会看到，式（5）₁中的$\nabla {\boldsymbol{v}_{\text{S}}}$不仅取决于骨架应变，还取决于固相基质应变和固相质量交换所产生的应变，因此式（5）₁中的$\nabla {\boldsymbol{v}_{\text{S}}}$与骨架应变并不是一一对应的关系。为了分析孔隙率对水-力-化学耦合的重要影响，需要利用式（2），（4）~（6）将固相应变分解为骨架应变、基质应变和固相质量交换所引起的应变。

利用Terzaghi有效应力$\tilde{\boldsymbol{\sigma}}=\boldsymbol{\sigma}-P_{\mathrm{F}} \boldsymbol{I}$和${\varphi _{\text{S}}} + {\varphi _{\text{F}}} = 1$，并结合式$\boldsymbol{\sigma}_{\mathrm{F}}=\varphi_{\mathrm{F}} P_{\mathrm{F}} \boldsymbol{I}$，$\boldsymbol{W}_{\mathrm{F}}=\boldsymbol{v}_{\mathrm{F}}-\boldsymbol{v}_{\mathrm{S}}$与式（2）可得

定义$\alpha$相材料的体应变为${\vartheta _\alpha } = \ln ({\rho _{{\text{R}}\alpha }}/{\rho _{{\text{R}}\alpha {\text{0}}}})$，${\rho _{{\text{R}}\alpha {\text{0}}}}$为$\alpha$相的初始材料密度，材料体应变${\vartheta _\alpha }$以受压为正。将${\vartheta _\alpha } = \ln ({\rho _{{\text{R}}\alpha }}/{\rho _{{\text{R}}\alpha {\text{0}}}})$代入式（7），则式（7）等式右侧第二项变为${P_{\text{F}}}{\varphi _{\text{S}}}{{d} ^{\text{S}}}{\rho _{{\text{RS}}}}/({\rho _{{\text{RS}}}}{d} t) =$${\varphi _{\text{S}}}{P_{\text{F}}}{{d} ^{\text{S}}}{\vartheta _{\text{S}}}/{d} t$，即此时固相材料（基质）体应变与流相孔压呈功共轭对，然而，固相材料（基质）体应变应该取决于固相基质压力。因此，需要将这一项的功共轭对进行变换以获得更符合多孔介质受力特征的能量平衡表达式。

将固相质量交换所引起的固相变形称为质量交换变形，记为${\varepsilon _{\text{c}}}$，其速率定义式为

定义${\boldsymbol{D}_{\text{H}}}$为

由式$\boldsymbol{D}_{\mathrm{S}}=\left[\nabla \boldsymbol{v}_{\mathrm{S}}+\left(\nabla \boldsymbol{v}_{\mathrm{S}}\right)^{\mathrm{T}}\right] / 2$，式(2)₁和（9）得

$\text{tr} {\boldsymbol{D}_{\text{H}}}$的物理含义是体应变速率，即式（10）表明与$\text{tr} {\boldsymbol{D}_{\text{H}}}$所对应的体积变形仅取决于固相体积分数，而固相体积分数与孔隙率一一对应，所以与${\boldsymbol{D}_{\text{H}}}$所对应的变形为骨架变形，${\boldsymbol{D}_{\text{H}}}$为固相骨架变形率。采用${\boldsymbol{D}_{\text{H}}}$来建立饱和多孔介质本构模型要比${\boldsymbol{D}_{\text{S}}}$更加简便。

令固相基质压力为${P_{\text{S}}} = ({P_{\text{T}}} - {\varphi _{\text{F}}}{P_{\text{F}}})/{\varphi _{\text{S}}}$，利用$\boldsymbol{D}_{\mathrm{S}}=\left[\nabla \boldsymbol{v}_{\mathrm{S}}+\left(\nabla \boldsymbol{v}_{\mathrm{S}}\right)^{\mathrm{T}}\right] / 2$，${\varphi _{\text{S}}} + {\varphi _{\text{F}}} = 1$，式（9）以及${\vartheta _\alpha } = \ln ({\rho _{{\text{R}}\alpha }}/{\rho _{{\text{R}}\alpha {\text{0}}}})$，则式（7）可表示为

式（11）表明，当采用${\boldsymbol{D}_{\text{H}}}$与有效应力$\tilde {\boldsymbol{\sigma}}$作为功共轭对时，与固相基质体应变${\vartheta _{\text{S}}}$呈共轭对的是固相基质压力，这与固相基质体应变取决于固相基质压力的力学规律一致，因而式（11）比式（7）更容易反映多孔介质各组分的力学规律。在试验中，固相基质压力无法直接测量，可通过测量总应力${P_{\text{T}}}$、孔压${P_{\text{F}}}$和孔隙率${\varphi _{\text{F}}}$并根据式${P_{\text{S}}}{\text{ = }}({P_{\text{T}}} - {\varphi _{\text{F}}}{P_{\text{F}}})/(1 - {\varphi _{\text{F}}})$来计算。

令$\boldsymbol{W}_\alpha=\boldsymbol{u}_\alpha-\boldsymbol{u}_{\mathrm{S}}$，将式（5）中的两个方程相加，再将式（6），（11）代入可得

其中$\hat{\boldsymbol{p}}_\alpha$为$\alpha$相的动量供给量。令: $\boldsymbol{F}_{\mathrm{S}}=\partial \boldsymbol{x} / \partial \boldsymbol{X}_{\mathrm{S}}$, $J_{\mathrm{S}}=\operatorname{det} \boldsymbol{F}_{\mathrm{S}}, \quad \boldsymbol{E}_{\mathrm{S}}=\left(\boldsymbol{F}_{\mathrm{S}}^{\mathrm{T}} \cdot \boldsymbol{F}_{\mathrm{S}}-\boldsymbol{I}\right) / 2, \quad J_{\mathrm{RS}}=\rho_{\mathrm{RS}} / \rho_{\mathrm{RS} 0}$, $J_{\mathrm{c}}=\exp \left(\int_0^t \hat{c}_{\mathrm{S}} / \rho_{\mathrm{S}} \mathrm{d} t\right), \quad \boldsymbol{F}_{\mathrm{H}}=\left[J_{\mathrm{c}} / J_{\mathrm{RS}}\right]^{1 / 3} \boldsymbol{F}_{\mathrm{S}}, \quad \boldsymbol{E}_{\mathrm{H}}=\left(\boldsymbol{F}_{\mathrm{H}}^{\mathrm{T}}\right.$. $\left.\boldsymbol{F}_{\mathrm{H}}-\boldsymbol{I}\right) / 2, \tilde{\boldsymbol{T}}_{\mathrm{H}}=\boldsymbol{F}_{\mathrm{H}}^{-1} \cdot \tilde{\boldsymbol{\sigma}} \cdot \boldsymbol{F}_{\mathrm{H}}^{-\mathrm{T}}$。其中$\boldsymbol{F}_{\mathrm{S}}$为固相变形梯度; $\boldsymbol{E}_{\mathrm{S}}$为固相Green应变; $\boldsymbol{F}_{\mathrm{H}}$为固相骨架的变形梯度; $\boldsymbol{E}_{\mathrm{H}}$为固相骨架Green应变; $\tilde{\boldsymbol{T}}_{\mathrm{H}}$为Piola-Kirchhoff有效应力。将$\tilde{\boldsymbol{T}}_{\mathrm{H}}$与$\boldsymbol{E}_{\mathrm{H}}$代入式(12) 可得

式（13）便是考虑化学活性时的盐溶液饱和多孔介质的能量平衡方程。

3.   熵流与熵产

根据式（9）可得

式（14）表明多孔固体总变形率${\boldsymbol{D}_{\text{S}}}$包括3个部分：固相基质体积变形率${{\text{d}}^{\text{S}}}{\vartheta _{\text{S}}}/{\text{d}}t$，固相骨架变形率${\boldsymbol{D}_{\text{H}}}$以及质量交换变形率${d} {\varepsilon _{\text{c}}}/{\text{d}}t$。式（13）表明${\boldsymbol{E}_{\text{H}}}$，${\vartheta _{\text{S}}}$与${\vartheta _{\text{F}}}$分别与${\tilde{ \boldsymbol{T}}_{\text{H}}}$，${\varphi _{\text{S}}}{P_{\text{S}}}$与${\varphi _{\text{F}}}{P_{\text{F}}}$呈功共轭关系。根据热力学理论，可选取${\boldsymbol{E}_{\text{H}}}$、${\vartheta _{\text{S}}}$和${\vartheta _{\text{F}}}$作为饱和多孔介质本构模型的应变状态变量；选取${{\tilde{\boldsymbol{T}}}_{\text{H}}}$、${\varphi _{\text{S}}}{P_{\text{S}}}$和${\varphi _{\text{F}}}{P_{\text{F}}}$作为应力状态变量。此外，在考虑化学活性时，还应选取各相各组分的质量分数及其化学势作为本构模型的化学状态变量。当本构模型需要描述塑性行为时，在内能表达式中还需选取塑性应变作为内变量，它的功共轭量称为耗散应力^[15]。在考虑化学活性的本构方程中，应选取两个塑性Green应变，分别为应力加载产生的骨架塑性应变$\boldsymbol{E}_{\text{H}}^{{\text{p1}}}$和化学加载产生的化学塑性应变$\boldsymbol{E}_{\text{H}}^{{\text{p2}}}$。假定固流两相具有相同温度，考虑到水溶剂或溶质与固相之间的吸附作用，流相的内能将受到骨架应变的影响。因此根据非平衡态热力学^[29]，式（13）中固相和流相的内能可分别表示为${\xi }_{\text{S}}（{\eta }_{\text{S}}, {\boldsymbol{E}}_{\text{H}}, {\boldsymbol{E}}_{\text{H}}^{\text{p1}}, {\boldsymbol{E}}_{\text{H}}^{\text{p2}}, {\vartheta }_{\text{S}})$和${\xi _{\text{F}}}({\eta _{\text{F}}}, {{\boldsymbol{E}}_{\text{H}}}, {\boldsymbol{E}}_{\text{H}}^{{\text{p1}}}, \boldsymbol{E}_{\text{H}}^{{\text{p2}}}, {\vartheta _{\text{F}}}, c_{\text{F}}^{\text{L}}, c_{\text{F}}^{\text{c}})$，其中${\eta _\alpha }$是$\alpha$相的熵密度，$c_{\text{F}}^{\text{L}} = \rho _{\text{F}}^{\text{L}}/{\rho _{\text{F}}}$与$c_{\text{F}}^{\text{c}} = \rho _{\text{F}}^{\text{c}}/{\rho _{\text{F}}}$分别为流相中溶剂与溶质的质量分数。由此可得内能的全微分为

根据热力学局部平衡假定，对比式（13）与（15）可得

令：

式中：$\bar \mu _{\text{F}}^{\text{L}}$，$\bar \mu _{\text{F}}^{\text{c}}$分别为与质量分数$c_{\text{F}}^{\text{L}}$与$c_{\text{F}}^{\text{c}}$呈功共轭关系的化学势；${\tilde{\boldsymbol{T}}}_{\text{H}}^{{\text{p1}}}$和${\tilde{\boldsymbol{T}}}_{\text{H}}^{{\text{p2}}}$为Piola-Kirchhoff耗散有效应力。根据试验成果获得内能表达式后，依据式（16a）~（16e）可获得有限应变条件下以内能形式表示的饱和多孔介质的弹性本构方程。

设$\alpha$相的Helmholtz自由能${\psi _\alpha } = {\xi _\alpha } - \theta {\eta _\alpha }$。令

将Helmholtz自由能定义式、式（13），（17）代入（15）并结合式（16a）~（16e）得

根据混合物理论^[17]有

令${\bar \mu _{\text{F}}} = c_{\text{F}}^{\text{L}}\bar \mu _{\text{F}}^{\text{L}} + c_{\text{F}}^{\text{c}}\bar \mu _{\text{F}}^{\text{c}}$，并将式（3）与（18）代入式（19），并利用$b\nabla \cdot {\boldsymbol{a}} = \nabla \cdot (b{\boldsymbol{a}}) - {\boldsymbol{a}} \cdot \nabla b$可得

根据热力学的熵分解理论^[29]，由式（20）可得熵流项为

熵产项为

在式（22）的右侧，第一项是热传导引起的熵产；第二项和第三项是塑性行为引起的熵产；第四项是渗流引起的熵产；第五项是由于在流体中扩散引起的熵产；第六项和第七项是由于固体和流体之间化学作用引起的熵产；最后一项是由于流相组分之间化学作用引起的熵产。根据热力学第二定律^[29]，熵产项应大于等于0。根据${{\boldsymbol{v}}_{\text{F}}} = (\rho _{\text{F}}^{\text{L}}{\boldsymbol{v}}_{\text{F}}^{\text{L}} + \rho _{\text{F}}^{\text{c}}{\boldsymbol{v}}_{\text{F}}^{\text{c}})/{\rho _{\text{F}}}$和${\boldsymbol{u}}_{\text{F}}^\beta = {\boldsymbol{v}}_{\text{F}}^\beta - {{\boldsymbol{v}}_{\text{F}}}$可得$\rho _{\text{F}}^{\text{L}}{\boldsymbol{u}}_{\text{F}}^{\text{L}} = - \rho _{\text{F}}^{\text{c}}{\boldsymbol{u}}_{\text{F}}^{\text{c}}$，将其代入式（22），即耗散率势函数为

4.   小应变条件下的一般本构方程

令$\boldsymbol{\varepsilon}_{\mathrm{S}}$与$\boldsymbol{\varepsilon}_{\mathrm{H}}$分别为固相应变张量$\boldsymbol{E}_{\mathrm{S}}$与骨架应变张量$\boldsymbol{E}_{\mathrm{H}}$在小应变条件下的近似值^[30]; $\boldsymbol{\varepsilon}_{\mathrm{H}}^{\mathrm{p} 1}$与$\boldsymbol{\varepsilon}_{\mathrm{H}}^{\mathrm{p} 2}$分别为骨架塑性应变张量$\boldsymbol{E}_{\mathrm{H}}^{\mathrm{p} 1}$与化学塑性应变张量$\boldsymbol{E}_{\mathrm{H}}^{\mathrm{p} 2}$在小应变条件下的近似值 ^[30]。小应变条件下有$\boldsymbol{D}_{\mathrm{S}}=\mathrm{d} \boldsymbol{\varepsilon}_{\mathrm{s}} / \mathrm{d} t$, $\boldsymbol{D}_{\mathrm{H}}=\mathrm{d} \boldsymbol{\varepsilon}_{\mathrm{H}} / \mathrm{d} t, \quad \rho_\alpha \approx \rho_{\alpha 0}, \quad \rho_{\mathrm{R} \alpha} \approx \rho_{\mathrm{R} \alpha 0}, \quad \tilde{\boldsymbol{T}}_{\mathrm{H}} \approx \tilde{\boldsymbol{\sigma}}$, $\tilde{\boldsymbol{T}}_{\mathrm{H}}^{\mathrm{p} 1} \approx \tilde{\boldsymbol{\sigma}}^{\mathrm{p} 1}$与$\tilde{\boldsymbol{T}}_{\mathrm{H}}^{\mathrm{p} 2} \approx \tilde{\boldsymbol{\sigma}}^{\mathrm{p} 2}$, 由式(14) 可得

式（24）表明，在小应变条件下，盐溶液饱和岩土材料的固相应变${\varepsilon _{\text{S}}}$可以分解为${{\boldsymbol{\varepsilon }}_{\text{H}}}$、${\vartheta _{\text{S}}}$以及${\varepsilon _{\text{c}}}$之和。

令${\varepsilon _{{\text{SV}}}}$与${\varepsilon _{{\text{HV}}}}$分别为固相体应变与固相骨架体应变，对式（24）取迹得

由式（10）和${\boldsymbol{D}_{\text{S}}} = \text{d} {{\boldsymbol{\varepsilon }}_{\text{S}}}/\text{d} t$可得在小应变条件下有${\varepsilon _{{\text{HV}}}} = \ln ({\varphi _{\text{S}}}/{\varphi _{{\text{S0}}}})$，注意到$\ln ({\varphi _{{\text{S0}}}}/{\varphi _{\text{S}}})$泰勒展开式中的高次项在小应变条件下可忽略，得

式中，${\varphi _{{\text{S0}}}}$为固相初始体积分数。同理，根据${\vartheta _\alpha } = \ln ({\rho _{{\text{R}}\alpha }}/{\rho _{{\text{R}}\alpha {\text{0}}}})$可得基质（材料）小应变表达式为

在小应变条件下, 式(16a)~(16e) 中的固相的内能变为$\xi_{\mathrm{S}}\left(\eta_{\mathrm{S}}, \boldsymbol{\varepsilon}_{\mathrm{H}}, \boldsymbol{\varepsilon}_{\mathrm{H}}^{\mathrm{p} 1}, \boldsymbol{\varepsilon}_{\mathrm{H}}^{\mathrm{p} 2}, \vartheta_{\mathrm{S}}\right)$, 流相的内能变为$\xi_{\mathrm{F}}\left(\eta_{\mathrm{F}}, \boldsymbol{\varepsilon}_{\mathrm{H}}, \boldsymbol{\varepsilon}_{\mathrm{H}}^{\mathrm{p} 1}, \boldsymbol{\varepsilon}_{\mathrm{H}}^{\mathrm{p} 2}, \vartheta_{\mathrm{F}}, c_{\mathrm{F}}^{\mathrm{L}}, c_{\mathrm{F}}^{\mathrm{c}}\right)$, 则以内能势函数表示的一般弹性本构方程式(16a)~(16e) 表示为

其中$\tilde{\boldsymbol{\sigma}}^{\mathrm{p} 1}$和$\tilde{\boldsymbol{\sigma}}^{\mathrm{p} 2}$是与$\boldsymbol{\varepsilon}_{\mathrm{H}}^{\mathrm{p} 1}$和$\boldsymbol{\varepsilon}_{\mathrm{H}}^{\mathrm{p} 2}$对应的耗散应力^{[15, 17]}。引入Helmholtz自由能为$\rho_{\alpha 0} \psi_\alpha=\rho_{\alpha 0} \xi_\alpha-\rho_{\alpha 0} \theta \eta_\alpha$, $\alpha \in\{\mathrm{S}, \mathrm{F}\}$, 此时固相Helmholtz自由能的表达式为$\psi_{\mathrm{S}}\left(\theta, \boldsymbol{\varepsilon}_{\mathrm{H}}, \boldsymbol{\varepsilon}_{\mathrm{H}}^{\mathrm{p} 1}, \boldsymbol{\varepsilon}_{\mathrm{H}}^{\mathrm{p} 2}, \vartheta_{\mathrm{S}}\right)$, 流相Helmholtz自由能表达式为$\psi_{\mathrm{F}}\left(\theta, \boldsymbol{\varepsilon}_{\mathrm{H}}, \boldsymbol{\varepsilon}_{\mathrm{H}}^{\mathrm{p} 1}, \boldsymbol{\varepsilon}_{\mathrm{H}}^{\mathrm{p} 2}, \vartheta_{\mathrm{F}}, c_{\mathrm{F}}^{\mathrm{L}}, c_{\mathrm{F}}^{\mathrm{c}}\right)$, 对$\rho_{\alpha 0} \psi_\alpha=\rho_{\alpha 0} \xi_\alpha-\rho_{\alpha 0} \theta \eta_\alpha$求全微分后将式(28) 代入得

由于固流两相的温度相同且恒定，自由能${\psi _{\text{S}}}$和${\psi _{\text{F}}}$表达式中的温度$\theta$可被省略。在大多数情况下，固相Helmholtz自由能可被简化为两部分：①固相材料体积变形所产生的自由能${\psi _{{\text{RS}}}}({\vartheta _{\text{S}}})$；②固相骨架变形所产生的自由能${\psi _{{\text{HS}}}}({\boldsymbol{\varepsilon}_{\text{H}}}, \boldsymbol{\varepsilon}_{\text{H}}^{{\text{p1}}}, \boldsymbol{\varepsilon}_{\text{H}}^{{\text{p2}}})$。根据式$c_{\text{F}}^{\text{L}} = 1 - c_{\text{F}}^{\text{c}}$，流相Helmholtz自由能${\psi _{\text{F}}}{\text{(}}{\boldsymbol{\varepsilon}_{\text{H}}}, \boldsymbol{\varepsilon}_{\text{H}}^{{\text{p1}}}, \boldsymbol{\varepsilon}_{\text{H}}^{{\text{p2}}}, {\vartheta _{\text{F}}}, c_{\text{F}}^{\text{L}}, c_{\text{F}}^{\text{c}})$可表达为${\psi _{\text{F}}}{\text{(}}{\boldsymbol{\varepsilon}_{\text{H}}}, \boldsymbol{\varepsilon}_{\text{H}}^{{\text{p1}}}, \boldsymbol{\varepsilon}_{\text{H}}^{{\text{p2}}}, {\vartheta _{\text{F}}}, c_{\text{F}}^{\text{c}})$，也可被简化为两部分：①流相基质（材料）的体积变形以及化学活性的作用而产生的自由能${\psi _{{\text{RF}}}}({\vartheta _{\text{F}}}, c_{\text{F}}^{\text{c}})$；②由于流相与固相之间的吸附作用而产生的自由能${\psi _{{\text{HF}}}}({\boldsymbol{\varepsilon}_{\text{H}}}, \boldsymbol{\varepsilon}_{\text{H}}^{{\text{p1}}}, \boldsymbol{\varepsilon}_{\text{H}}^{{\text{p2}}})$。令${\varPsi }_{\text{HF}}(\boldsymbol{\varepsilon}_{\text{H}}, \boldsymbol{\varepsilon}_{\text{H}}^{\text{p1}}, \boldsymbol{\varepsilon}_{\text{H}}^{\text{p2}})={\rho }_{\text{F0}}{\psi }_{\text{HF}}(\boldsymbol{\varepsilon}_{\text{H}}, \boldsymbol{\varepsilon}_{\text{H}}^{\text{p1}}, \boldsymbol{\varepsilon}_{\text{H}}^{\text{p2}})，{\varPsi }_{\text{HS}}(\boldsymbol{\varepsilon}_{\text{H}}, \boldsymbol{\varepsilon}_{\text{H}}^{\text{p1}}, \boldsymbol{\varepsilon}_{\text{H}}^{\text{p2}})=$${\rho }_{\text{S0}}{\psi }_{\text{HS}}(\boldsymbol{\varepsilon}_{\text{H}}, \boldsymbol{\varepsilon}_{\text{H}}^{\text{p1}}, \boldsymbol{\varepsilon}_{\text{H}}^{\text{p2}})，{\varPsi }_{\text{RS}}({\vartheta }_{\text{S}})={\rho }_{\text{RS0}}{\psi }_{\text{RS}}({\vartheta }_{\text{S}})，{\varPsi }_{\text{RF}}({\vartheta }_{\text{F}},$$c_{\text{F}}^{\text{c}}) = {\rho _{{\text{RF0}}}}{\psi _{{\text{RF}}}}({\vartheta _{\text{F}}}, c_{\text{F}}^{\text{c}})$，根据上述分析由式（29）可得

式（30）是在小应变条件下以自由能势函数表示的一般弹性本构方程。由于考虑到温度恒定，则耗散率势函数式（23）在小应变条件下可表示为

当以下关系成立时，式（31）显然成立。

Collins等^[20]认为，一般可假定耗散势存在：

式（33a）~（33e）基于热力学第二定律的耗散势建立，反映了多孔介质非弹性、扩散和化学反应等不可逆过程，其塑性本构方程与渗流-扩散控制方程可据此获得。

5.   算例验证

上面所建立的是盐溶液饱和岩土的一般本构规律。具体建模时，在上述一般本构规律指导下，可以结合黏土的微细观结构，从物理意义的角度出发确定具体材料的本构方程；也可以结合试验成果，通过拟合试验数据来建立具体材料的本构方程。本文结合颜荣涛等^[5]开展的NaCl溶液饱和黏土的等向压缩试验来建立饱和黏土的等向压缩本构模型。

为了验证上述盐溶液饱和岩土的本构框架的实用性，现结合具体试验成果来建立NaCl溶液饱和黏土的本构建模。试验用土取自南宁的五象新区某工地，为高液限黏土^[5]。黏土试样由NaCl溶液饱和，试验过程为等向加压到1200 kPa后卸压回弹至100 kPa。结合试样在去离子水和质量分数为3%，11%的NaCl溶液饱和下的等向压缩试验数据，确定模型参数。利用该模型参数模拟了试样在质量分数分别为0.6%和6%的NaCl溶液饱和时的等向压缩力学行为。

在等向压缩试验中，剪应力为0，试样在平均有效应力作用下只产生体应变，因此式（30）可退化为

式中：$\tilde P$为平均有效应力；${\tilde P^{{\text{p1}}}}$和${\tilde P^{{\text{p2}}}}$为反映第1和第2屈服面的耗散有效应力；$\varepsilon _{{\text{HV}}}^{{\text{p1}}}, \varepsilon _{{\text{HV}}}^{{\text{p2}}}$分别为骨架塑性体应变与化学塑性体应变。式（34）表明对于饱和多孔介质的本构关系包含固相和流相本构关系，这与Biot固结理论以及陈正汉的研究成果相一致^{[18-19, 28]}。同时式（34）中与化学活性相关的状态变量是无量纲的质量分数而不是摩尔浓度，这与当前采用摩尔浓度作为状态变量的盐溶液饱和岩土本构模型有所区别。此外，化学势在试验中不易定量，因此采用质量分数作为化学变量来建立土体固相本构模型则更为方便。

根据以往研究可知，对于大多数弹塑性问题，耗散有效应力${\tilde P^{{\text{p1}}}}$和${\tilde P^{{\text{p2}}}}$与Terzaghi有效应力大小相等${\tilde P^{{\text{p1}}}} = {\tilde P^{{\text{p2}}}} = \tilde P$，因此${\tilde P^{{\text{p1}}}}$和${\tilde P^{{\text{p2}}}}$也被称为Terzaghi耗散有效应力^[15-16]。

目前在关于盐溶液饱和黏土的实际工程应用与研究中，通常不考虑固相与液相之间以及流相中各组分之间的质量交换。因此有${\hat c_{F} } = {\hat c_{\text{S}}} = 0$，$\hat c_{\text{F}}^{\text{L}} = \hat c_{\text{F}}^{\text{c}} = 0$与${\varepsilon _{\text{c}}} = 0$。则根据式（25）可得土体固相的体应变为

在等向压缩试验中，（33a）~（33c）可写为

式（36）表明盐溶液饱和黏土的本构关系还应包括塑性本构关系与流体的渗流规律以及溶质的扩散方程。

5.1   固体基质（土颗粒）的本构关系

在土力学中，固相基质的变形即为土体颗粒的变形。当土骨架破坏时，一般认为固相基质仅产生弹性应变。对于线弹性模型，固相基质的Helmholtz自由能由式（34）₄可表示为

式中，${K_{{\text{RS}}}}$为土颗粒的体积模量，将式（37）代入式（34）₄得：

5.2   土骨架的弹性本构关系

设骨架弹性体应变为$\varepsilon _{{\text{HV}}}^{\text{e}} = {\varepsilon _{{\text{HV}}}} - \varepsilon _{{\text{HV}}}^{{\text{p1}}} - \varepsilon _{{\text{HV}}}^{{\text{p2}}}$，根据本文理论框架，重新整理试验数据可以得到土体骨架弹性应变^[5]，如图 1所示。土骨架的Helmholtz自由能${\varPsi _{{\text{HS}}}}({\varepsilon _{{\text{HV}}}}, \varepsilon _{{\text{HV}}}^{{\text{p1}}}, \varepsilon _{{\text{HV}}}^{{\text{p2}}})$与${\varPsi _{{\text{HF}}}}({\varepsilon _{{\text{HV}}}}, \varepsilon _{{\text{HV}}}^{{\text{p1}}}, \varepsilon _{{\text{HV}}}^{{\text{p2}}})$可表示为$\varepsilon _{{\text{HV}}}^{\text{e}}$的函数^[16]，结合式（34），进一步可表示为

图  1  骨架弹性体应变随平均有效应力变化

Figure  1.  Variation of elastic volumetric strain of skeleton with average effective stress

下载: 全尺寸图片 幻灯片

式中：${\tilde P_{\text{0}}}$为初始有效应力；${\lambda ^{\text{e}}}$为回弹曲线斜率；试验结果表明，随着NaCl溶液的质量分数增加，${\lambda ^{\text{e}}}$基本保持不变^[5]。${\varepsilon _{{\text{Hmax}}}}$为试样去离子水饱和条件下，平均有效应力为${\tilde P_{\text{0}}}$时由于吸附作用而产生的最大膨胀应变。饱和土的体积分数与含水率w之间的关系为${\varphi _{\text{S}}} = {d_{\text{F}}}/(1 + w{d_{\text{S}}})$，${d_{\text{F}}}$与${d_{\text{S}}}$分别是流相和固相的相对密度，因此最大膨胀应变${\varepsilon _{{\text{Hmax}}}}$与最大含水率${w_{\max }}$之间关系为${\varepsilon _{{\text{Hmax}}}} = ({\varphi _{{\text{S0}}}} + {\varphi _{{\text{S0}}}}{d_{\text{S}}}{w_{\max }} - {d_{\text{F}}})/({\varphi _{{\text{S0}}}} + {\varphi _{{\text{S0}}}}{d_{\text{S}}}{w_{\max }})$。将式（39）代入式（34）₁可得

式（40）还可进一步表示为

5.3   土骨架的塑性本构关系

根据试验结果^[5]和骨架的弹性体应变，可以得到骨架塑性体应变随着平均有效应力的变化，如图 2所示。根据图 2，并结合式（36）₁，${\Theta _{{\text{p1}}}}$可表示为

图  2  骨架塑性体应变随平均有效应力变化

Figure  2.  Variation of plastic volumetric strain of skeleton with average effective stress

下载: 全尺寸图片 幻灯片

将式（42）代入式（36）₁得

为了使（32a）₁成立，即当$\text{d} \varepsilon _{{\text{HV}}}^{{\text{p1}}}/\text{d} t \geqslant 0$时，$\text{d} \tilde P =$ $\text{d} {\tilde P^{{\text{p1}}}} > 0$，式（43）可进一步表示为

式中：$\varepsilon _{{\text{HV0}}}^{{\text{p1}}}(c_{\text{F}}^{\text{c}})$为质量分数为$c_{\text{F}}^{\text{c}}$的NaCl溶液饱和黏土骨架的初始塑性体应变；${\lambda ^{\text{p}}}(c_{\text{F}}^{\text{c}})$是塑性压缩曲线斜率，它与溶液质量分数有关，如图 3所示。塑性压缩曲线斜率与溶液质量分数的关系为

图  3  塑性压缩线斜率随NaCl质量分数变化

Figure  3.  Variation of slope of plastic compression line with NaCl mass fraction

下载: 全尺寸图片 幻灯片

式中，$\lambda _1^{\text{p}}$，$\lambda _2^{\text{p}}$以及$\lambda _3^{\text{p}}$为塑性压缩参数，可根据试验结果拟合获得。与式（45）不同，Loret等^[13]采用化学势$\bar \mu _{\text{F}}^{\text{c}}$代替质量分数$c_{\text{F}}^{\text{c}}$来建立塑性压缩参数公式。

图 4为在$\varepsilon _{{\text{HV}}}^{{\text{p1}}}{\text{ - }}\ln \tilde P$平面上盐溶液饱和黏土的塑性压缩曲线示意图，类比于孙德安等^[31]开展的吸力对非饱和土压缩曲线影响的研究，可认为不同质量分数的NaCl溶液饱和黏土的塑性压缩曲线相交于一点。其中${\tilde P_{\text{c}}}(c_{\text{F}}^{\text{c}})$和${\tilde P_{\text{c}}}(0)$分别为NaCl溶液饱和与去离子水饱和黏土的屈服应力；${\tilde P_{{\text{ref}}}}$与$\varepsilon _{{\text{Href}}}^{{\text{p1}}}$分别为参考应力与应变；当初始有效应力${\tilde P_{\text{0}}}{\text{ = }}{\tilde P_{{\text{ref}}}}$，则有$\varepsilon _{{\text{HV0}}}^{{\text{p1}}}(c_{\text{F}}^{\text{c}}) = \varepsilon _{{\text{Href}}}^{{\text{p1}}}$。根据式（44）和图 4可得

图  4  在$\varepsilon _{{\text{HV}}}^{{\text{p1}}}{\text{ - }}\ln \tilde P$平面上盐溶液饱和黏土的压缩曲线

Figure  4.  Compression curves of clay saturated with salt solution in $\varepsilon _{{\text{HV}}}^{{\text{p1}}}{\text{ - }}\ln \tilde P$ plane

下载: 全尺寸图片 幻灯片

即可得屈服应力为

将式（45）代入式（47），得到的计算结果与试验数据进行对比，如图 5所示。图 5表明，实测数据与计算结果基本一致。

图  5  屈服应力的模拟结果与试验数据对比

Figure  5.  Comparison of simulated results and experimental data for yield effective stress

下载: 全尺寸图片 幻灯片

从饱和黏土的等向压缩试验数据可以看出，在相同的等向应力荷载作用下，溶液质量分数的变化（即化学加载）会产生化学塑性体应变$\varepsilon _{{\text{HV}}}^{{\text{p2}}}(c_{\text{F}}^{\text{c}})$，此塑性变形可采用第二个屈服面来模拟。根据NaCl溶液饱和黏土等向压缩试验所体现的化学塑性变形变化规律^[5]，第二屈服面的耗散势函数表示为

其中，$\varepsilon _{\text{p}}^{\text{c}}$与${\lambda ^{\text{c}}}$为塑性应变参数，把式（48）代入式（36）₂得

根据式（49）可得

在等向压缩试验中，${\tilde P^{{\text{p2}}}} = \tilde P > 0$，由式（32a）₂所要求的不等式${\tilde P^{{\text{p2}}}}\text{d} \varepsilon _{{\text{HV}}}^{{\text{p2}}} \geqslant 0$可得第二屈服面的塑性变形存在的条件是$\varepsilon _{{\text{HV}}}^{{\text{p2}}}$需满足${\text{d}}\varepsilon _{{\text{HV}}}^{{\text{p2}}} \geqslant 0$。基于等向压缩试验成果^[5]，可以获得化学塑性体应变随NaCl溶液质量分数的变化规律，如图 6所示。

图  6  化学塑性体应变随NaCl溶液质量分数变化

Figure  6.  Variation of chemical plastic strain with mass fraction of NaCl solution

下载: 全尺寸图片 幻灯片

5.4   土体固相的本构模型

根据式（35），（38），（42），（44），（50）可得NaCl溶液饱和黏土的固相应变为

参数${\lambda ^{\text{c}}}$与$\lambda$取值情况如表 1所示，其中$c_{{\text{Fc}}}^{\text{c}}$为化学加载达到塑性屈服的临界值。

表  1  参数${\lambda ^{\text{c}}}$与$\lambda$取值情况

Table  1.  Values of parameters${\lambda ^{\text{c}}}$and$\lambda$

取值	$c_{\text{F}}^{\text{c}} < c_{{\text{Fc}}}^{\text{c}}$ $\tilde P < {\tilde P_{\text{c}}}$	$c_{\text{F}}^{\text{c}} \geqslant c_{{\text{Fc}}}^{\text{c}}$ $\tilde P < {\tilde P_{\text{c}}}$	$c_{\text{F}}^{\text{c}} \geqslant c_{{\text{Fc}}}^{\text{c}}$ $\tilde P \geqslant {\tilde P_{\text{c}}}$	$c_{\text{F}}^{\text{c}} < c_{{\text{Fc}}}^{\text{c}}$ $\tilde P \geqslant {\tilde P_{\text{c}}}$
$\lambda$	${\lambda ^{\text{e}}}$	${\lambda ^{\text{e}}}$	${\lambda ^{\text{e}}} + {\lambda ^{\text{p}}}(c_{\text{F}}^{\text{c}})$	${\lambda ^{\text{e}}} + {\lambda ^{\text{p}}}(c_{\text{F}}^{\text{c}})$
${\lambda ^{\text{c}}}$	0	${\lambda ^{\text{c}}}$	${\lambda ^{\text{c}}}$	0

下载: 导出CSV 

| 显示表格

利用去离子水、3%与11% NaCl溶液饱和黏土试样的试验数据来标定模型参数，由此获取的模型参数如表 2所示。利用表 2的模型参数和本构方程理论模拟0.6%和6%试样的应变特性，并与试验数据相对比以验证本文本构方程的合理性，模拟结果与试验数据对比如图 7所示。从图 7可以看出该模型能有效地反映NaCl溶液对黏性土等向压缩行为的影响，从而表明采用混合物理论^[17]与经典热力学溶液理论^[26]所建立的盐溶液饱和岩土的本构理论可以指导盐溶液饱和黏土固相本构模型的建立。

表  2  等向压缩试验的模型参数选取

Table  2.  Model parameters for isotropic compression tests

$\lambda _1^{\text{p}}$	$\lambda _2^{\text{p}}$	$\lambda _3^{\text{p}}$	${\lambda ^{\text{c}}}$	${\lambda ^{\text{e}}}$
0.0105	-0.147	0.0099	1/0.09	0.014
$\varepsilon _{\text{p}}^{\text{c}}$	${\varepsilon _{{\text{Hmax}}}}$	${\tilde P_{{\text{ref}}}}$/kPa	${\tilde P_{\text{c}}}(0)$/kPa	${K_{{\text{RS}}}}$/GPa
-0.08	0.02	5	75	20

下载: 导出CSV 

| 显示表格

图  7  等向固结试验数据与模拟结果对比

Figure  7.  Comparison of experimental data of isotropic consolidation with simulation results

下载: 全尺寸图片 幻灯片

5.5   流相材料的应力应变关系

在以往研究中，流相基质本构关系也称为流相材料的本构关系^[15-16]。本文假定流相为无黏性流体，流相材料Gibbs自由能为${G_{{\text{RF}}}}({\vartheta _{\text{F}}}, c_{\text{F}}^{\text{c}}) = {\varPsi _{{\text{RF}}}}({\vartheta _{\text{F}}}, c_{\text{F}}^{\text{c}}) - {P_{\text{F}}}{\vartheta _{\text{F}}}$，根据式（34）₅和式（34）₆得

根据式（52）可得

根据经典热力学溶液理论^[26]可知，溶质与溶剂的化学势之差可表示为

式中：$\bar \mu _{\text{F}}^{{\text{c0}}}(\theta , {P_0})$与$\bar \mu _{\text{F}}^{{\text{L0}}}(\theta , {P_0})$分别为标准状态下NaCl溶质与水溶剂的化学势；$M_{\text{F}}^{\text{c}}$与$M_{\text{F}}^{\text{L}}$分别为NaCl溶质与水溶剂的的摩尔质量；R为气体通用常数。通常认为理想溶液的偏摩尔容积$- {\bar V_{\text{F}}}$保持不变^[26]，对式（54）关于${P_{\text{F}}}$求偏导，将其代入式（53）可得

式中，${K_{\text{F}}}$为NaCl溶液的体积模量，图 8给出了NaCl溶液的体应变随着NaCl质量分数的变化关系^[32]。

图  8  NaCl溶液体应变随NaCl质量分数变化

Figure  8.  Variation of strain of NaCl solution with NaCl mass fraction

下载: 全尺寸图片 幻灯片

根据式（55）可得

式（56）表明溶液体应变${\vartheta _{\text{F}}}$是由孔隙水压力${P_{\text{F}}}$与NaCl质量分数$c_{\text{F}}^{\text{c}}$所共同决定，其中任何一个因素改变都将会影响流相材料的应变。

5.6   流相渗流的达西定律与一维固结控制方程

在等向压缩试验中，试样在上下两端双面排水，因此可认为流相主要发生竖直方向的渗流。令式（37）₂ ${\Theta _{\text{w}}}({W_{\text{F}}}){\text{ = }}{\gamma _{\text{w}}}{\varphi _{{\text{F0}}}}{\varphi _{\text{F}}}W_{\text{F}}^2/(2{k_{\text{F}}})$，其中，${k_{\text{F}}}$为渗透系数，${\gamma _{\text{w}}}$为水的重度。将${\Theta _{\text{w}}}({W_{\text{F}}})$代入式（37）₂可得

NaCl溶液的动量平衡方程为^[10]

将式（58）代入式（59）可得

在小应变条件下，由式（17）有

根据$\partial ({\rho _{{\text{F0}}}}{\xi _{\text{F}}})/\partial {\varepsilon _{{\text{HV}}}} = \partial {\varPsi _{{\text{HF}}}}/\partial {\varepsilon _{{\text{HV}}}}$，$\partial ({\rho _{{\text{F0}}}}{\xi _{\text{F}}})/$$\partial \varepsilon _{{\text{HV}}}^{{\text{p1}}}$$= \partial {\varPsi _{{\text{HF}}}}/\partial \varepsilon _{{\text{HV}}}^{{\text{p1}}}$，$\partial ({\rho _{{\text{F0}}}}{\xi _{\text{F}}})/\partial \varepsilon _{{\text{HV}}}^{{\text{p2}}} = \partial {\varPsi _{{\text{HF}}}}/\partial \varepsilon _{{\text{HV}}}^{{\text{p2}}}$，再结合$\varepsilon _{{\text{HV}}}^{\text{e}} = {\varepsilon _{{\text{HV}}}} - \varepsilon _{{\text{HV}}}^{{\text{p1}}} - \varepsilon _{{\text{HV}}}^{{\text{p2}}}$，式（60）可表示为

根据式（39）~（41）以及（50）得

将式（61），（62）代入式（59）可得

式（63）为广义达西定律。当不考虑由于溶质与水的吸附作用而产生的膨胀变形时，有${\varepsilon _{{\text{Hmax}}}} + \varepsilon _{{\text{HV}}}^{{\text{p2}}}(c_{\text{F}}^{\text{c}}) = 0$，式（63）即可退化为经典的达西定律。

在实际工程中，岩土工程界关心的是流体从多孔固体中的流出或流入量。根据混合物理论和文献[16]，可得流体渗入量为

将${\varphi _{{\text{F0}}}} + {\varphi _{{\text{S0}}}} = 1$与式（38），（51）与（55）代入式（64）可得

利用${\zeta _{\text{F}}} = {\varphi _{{\text{F0}}}}({\varepsilon _{{\text{SV}}}} - {\varepsilon _{{\text{FV}}}})$和${W_{\text{F}}} = {v_{\text{F}}} - {v_{\text{S}}}$可得$\partial ({\varphi _{{\text{F0}}}}{W_{\text{F}}})/$ $\partial z = - \partial {\zeta _{\text{F}}}/\partial t$，再结合式（63）和（65）可得

式（66）是考虑固体基质（土颗粒）可压缩和化学活性的饱和黏土一维固结控制方程。

5.7   NaCl溶质的渗流-扩散方程

令${\Theta _{\text{d}}}(\rho _{\text{F}}^{\text{c}}u_{\text{F}}^{\text{c}}) = {k_{\text{d}}}{(\rho _{\text{F}}^{\text{c}}u_{\text{F}}^{\text{c}})^2}/2$，式（37）₃可表示为

其中${k_{\text{d}}}$为扩散系数，z为扩散距离。在小应变条件下，有${v_{\text{S}}} \approx 0$，${\rho _{\text{F}}} \approx {\rho _{{\text{F0}}}}$，将式（67）代入式（3）₂，并利用式$v_{\text{F}}^{\text{c}} = u_{\text{F}}^{\text{c}} + {v_{\text{F}}} = u_{\text{F}}^{\text{c}} + {W_{\text{F}}} + {v_{\text{S}}}$可得

将式（54）代入式（68）可得

式（69）即为NaCl溶质在饱和黏土中的渗流-扩散方程。

5.8   渗透压

为了使式（69）在${W_{\text{F}}} = 0$时$\partial c_{\text{F}}^{\text{c}}/\partial t = 0$成立，即渗流量为零时质量分数保持不变，需满足以下关系：

根据热力学溶液理论^[26]，对于理想的NaCl溶液，其溶剂的化学势可表示为

式中，${\bar V_{{\text{m, L}}}}$为水溶剂的偏摩尔容积。将式（54）代入式（71）得

半透膜阻止溶质穿过半透膜扩散，相当于NaCl溶质在半透膜上承受一个阻力使其达到化学势平衡，即$b_{\text{F}}^{\text{c}} = - \partial \bar \mu _{\text{F}}^{\text{c}}/\partial z$。当外载荷满足$b_{\text{F}}^{\text{c}} = - \partial \bar \mu _{\text{F}}^{\text{c}}/\partial z$，则相当于半透膜作用。将式$b_{\text{F}}^{\text{c}} = - \partial \bar \mu _{\text{F}}^{\text{c}}/\partial z$与式（72）代入式（70）得

对式（73）积分，先利用${P_{\text{F}}} + \int_0^z {(b_{\text{F}}^{\text{L}}M_{\text{F}}^{\text{L}}/} {\bar V_{{\text{m, L}}}})$ $\text{d} z = 0$，再根据$c_{\text{F}}^{\text{c}}$非常小时可略去对数函数泰勒展开的高次项，最后采用$c_{\text{F}}^{\text{c}} = \rho _{\text{F}}^{\text{c}}/{\rho _{\text{F}}}$代替$c_{\text{F}}^{\text{c}}$，可得水溶剂化学势平衡的表达式为

式中：$\rho _{\text{F}}^c/M_{\text{F}}^{\text{c}}$，$\rho _{\text{F}}^{\text{L}}/M_{\text{F}}^{\text{L}}$分别为溶质与溶剂的摩尔数，故式（74）即为广义渗透压的一般表达式。广义渗透压不仅与溶质的质量分数有关，还与所受外荷载有关。当流相动量达到平衡${W_{\text{F}}} = 0$和流相中组分的化学势达到热力学平衡$\partial c_{\text{F}}^{\text{c}}/\partial t = 0$时，式（74）就变成了经典渗透压公式。饱和土在外荷载作用下会产生固结，此时${W_{\text{F}}} \ne 0$，故经典渗透压公式在饱和多孔介质固结过程中不成立。

6.   结论

（1）把固相变形分解为骨架变形、固相基质变形和质量交换变形，以凸显孔隙率和化学活性在水-力-化学多场耦合中的关键作用；采用现时质量分数作为化学状态变量来反映饱和岩土孔隙溶液的化学活性，改变了以往采用渗透压和溶质摩尔浓度作为化学状态变量建立模型的传统思路。

（2）研究表明，弹性本构方程可以通过自由能势函数推导而得；塑性本构方程、广义达西定律可以通过耗散势函数推导而得。采用混合物理论和热力学溶液理论建立的盐溶液饱和多孔介质本构理论框架可以用来指导创建本构模型。

（3）在上述理论指导下，结合NaCl溶液饱和黏土的等向压缩试验成果，建立了考虑水-力-化学耦合特性的NaCl溶液饱和黏土本构模型。该本构模型包括固流两相本构方程和溶质渗流-扩散方程。经典渗透压公式需同时满足流相动量守恒方程和热力学平衡条件才成立。