Thermodynamic model of coupled temperature and pressure effects for hydrate-bearing sediments within particle rearrangement theory
-
摘要: 水合物沉积物赋存时的温度和孔隙水压变化对沉积物的力学性质影响明显。基于颗粒重组理论，通过引入温-压耦合系数δ、剪应变和体应变胶结应力衰减方式(λv，λs)和水合物饱和度影响下的剪胀方程d，建立了水合物沉积物温-压耦合热力学模型。利用数值模拟结果同室内试验结果的对比分析, 从宏观力学性质变化和微观作用机理两个角度探讨了围压、水合物饱和度和温-压耦合系数对沉积物力学特性的影响。最后，对模型中的刚度系数和体应变衰减参数进行敏感性分析。结果表明：引入温-压耦合系数的水合物沉积物模型能够很好的描述沉积物赋存的温度和水压同其力学特性变化关系。环境温度降低和水压增大增强了微观层面水合物的胶结强度和刚度，宏观表现为峰值强度、应变软化和剪胀性提高。刚度系数γ通过增大沉积物初始刚度使得沉积物表现出更强的峰值强度。体应变衰减参数λv通过提高胶结应力衰减速率增强沉积物的应变软化特性。Abstract: The mechanical properties of hydrate-bearing sediments are significantly influenced by changes in temperature and water pressure. Based on the particle rearrangement theory, a thermodynamic model that couples temperature and pressure, incorporating the dilation equation considering the effects of hydrate saturation and introducing the bond degradation parameters caused by shear and volumetric strains, is developed to describe the mechanical behaviors of hydrate-bearing sediments. The numerical simulation results are compared with the laboratory experiments to explore the effects of confining pressure, hydrate saturation and the temperature-pressure coupling coefficient on the mechanical properties of sediments from both macroscopic and microscopic perspectives. Finally, the sensitivity analyses are conducted on the stiffness coefficient and bond degradation parameter. The results indicate that the introduction of a temperature-pressure coupling coefficient in the model effectively describes the relationship between the mechanical properties of sediments and temperature and water pressure during deposition. Decreasing environmental temperature and increasing water pressure enhance the bond strength and stiffness of hydrates at the microscopic level, resulting in the increased peak strength, strain softening and shear dilation at the macroscopic level. Increasing the stiffness coefficient γ enhances the peak strength of sediments by increasing the initial stiffness of sediments. The bond degradation parameter enhances the strain-softening behaviors of sediments by increasing the rate of bond degradation.
-
0. 引言
水合物是由甲烷气体和液态水在低温和高压条件下形成的固态笼型晶体物质，具有能量密度高、分布广、储藏大和燃烧清洁的优点[1-2]。水合物常以胶结、填充、包裹和支撑的形式存在海底沉积层和冻土区域土颗粒之间的孔隙中，对沉积层土体强度、刚度和变形特性影响显著[3]。当前，对水合物的研究逐步从资源勘探过渡到勘查试开采阶段[4]。开采方式包括升温法、降压法、置换法和联合法，其目的是使得水合物所处的环境温度和水压力条件逼近对应的相平衡点，从而使水合物发生分解，实现水合物的开采[5]。然而，在水合物试开采阶段，沉积层中水合物的分解会引起沉积层温度、水压力、颗粒之间的接触应力和胶结应力在空间和时间尺度分布不均衡，对开采平台和海底环境产生较大威胁[6-7]。另一方面，学者们对含水合物沉积物的力学性质研究集中在水合物饱和度、围压和赋存模式等影响因素上，对于温度和孔隙压力的影响研究较少[8-9]。因此，有必要对沉积物赋存时的温度和水压力变化引起的力学性质演化进行研究。
近年来，学者们通过人工合成水合物沉积物试样, 开展室内试验研究水合物赋存条件对沉积物力学性质的影响[10-11]。研究表明在水合物未发生相变区域，水合物沉积物的宏观力学性质表现为赋存环境温度越低，孔隙压力和有效围压越大, 其对应的强度越高，应变软化性能和剪胀特性越明显[11]。而在微观性质研究上，Winters等[12]利用声波检测装置对砂土试样进行三轴剪切时发现孔隙中的游离气体对试样的强度响应有延迟作用。Chaouachi等[13]提出在微观层次, 相邻土颗粒之间的连接存在胶结模式和吸附水膜两种模式。赋存环境的变化同时影响着水合物产生量和吸附水膜的厚度。吸附水膜越薄, 离子强度越高，界面静电势越低，吸附水膜与土颗粒和水合物之间的范德华力增强[14]。宏观表现为沉积物试样强度增大。蒋明镜等[15]提出在温度和孔隙压力坐标系内，以状态点至相变线的最近距离为温-压耦合参数。该理论能够综合的考虑试验时的温度和孔隙压力对土体强度的耦合影响，但是存在着数值计算复杂问题。另一方面，文献[16]提出在温度和孔隙压力坐标内，采用试验温度和对应孔隙压力的相平衡温度差作为温-压耦合参数引入模型计算中。然而该方法忽略了试验孔隙水压对沉积物强度的影响, 对于不同孔隙压力和不同温度条件下的试样强度预测存在困难。因此，有必要定义一个同时考虑温度和水压影响且计算方式简单，物理含义明确的温-压耦合系数。
随着对水合物沉积物的力学特性的深入了解，研究者们开始采用不同的理论方法建立水合物沉积物的本构关系[17-18]。其中基于经典弹塑性理论的沉积物本构关系应用最为广泛[19-20]。然而，这类模型在建立过程中，对于流动法则和屈服函数的建立基于经验公式，缺少严格的理论论证。近年来，学者们开始采用颗粒热动力学理论建立材料的本构关系。该方法是基于颗粒重组和混合物理论，采用迁移系数矩阵和能量密度函数表征颗粒类材料能量耗散过程，能够建立材料的宏观力学行为同微观特性之间的联系，在非饱和土和特殊土的本构关系构建中得到成功应用[21-22]。
以水合物沉积物的颗粒重组理论为基础建立了温-压耦合作用下水合物沉积物的热力学模型。该模型定义了温-压耦合系数表达式，并将其引入到胶结应力和刚度参数表达式中。在理论模型中考虑水合物饱和度影响下的剪胀方程。通过对水合物饱和度、围压和温-压耦合系数影响下沉积物的力学特性进行预测，并将模拟结果与试验结果进行对比分析来验证模型的有效性。最后，对建模过程中的刚度参数和体应变衰减参数进行敏感性分析，从颗粒重组角度对参数的物理含义进行确定。
1. 热力学理论模型建立
颗粒重组理论是由Jiang等[23]研究单相颗粒物质材料的本构关系时提出的。该模型首次引入颗粒熵和颗粒温度的概念，并建立了颗粒材料宏观力学行为和细观行为的联系，是一种多尺度模型。在此基础上，通过引入土水特征曲线并考虑饱和度变化引起的能量耗散，Yang等[21]利用该理论建立了非饱和土的热-水-力耦合模型。模型能够反映非饱和土的孔隙压力变化、应变硬化（软化）和剪胀特性。近年来，Bai等[22]对颗粒热动力学模型参数的物理含义、量纲和取值范围进行了标定，为特殊土的本构关系建立提供理论依据，使得颗粒热动力学模型趋于成熟。颗粒热动力学理论的经典表达式为
$$ \pi_{i j}=\frac{\partial\left[W_{\mathrm{e}}\left(B, \varepsilon_{\mathrm{v}}^{\mathrm{e}}, \varepsilon_{\mathrm{s}}^{\mathrm{e}}, T, C_{\mathrm{s}}\right)\right]}{\partial \varepsilon_{i j}^{\mathrm{e}}} , $$ (1) $$ d_{\mathrm{t}} \varepsilon_{i j}^{\mathrm{p}}=C_1\left(T_{\mathrm{GG}}\right)^{0.5}\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}+\frac{C_{\mathrm{s}}}{K_{\mathrm{e}}}\right) \delta_{i j}+\left(T_{\mathrm{GG}}\right)^{0.5} e_{i j}^{\mathrm{e}} , $$ (2) $$ d_{\mathrm{t}} T_{\mathrm{GG}}=C_2 C_5 \frac{d_{\mathrm{t}} \varepsilon_{\mathrm{s}} d_{\mathrm{t}} \varepsilon_{\mathrm{s}}}{\rho^{\mathrm{s}}}+C_3 C_5 \frac{d_{\mathrm{t}} \varepsilon_{\mathrm{v}} d_{\mathrm{t}} \varepsilon_{\mathrm{v}}}{\rho^{\mathrm{s}}}+ $$ $$ C_4 \frac{R d_{\mathrm{t}} \varphi}{\rho^{\mathrm{s}}}-C_4 \frac{T_{\mathrm{GG}}}{\rho^{\mathrm{s}}}。 $$ (3) 式中：πij为土颗粒之间的接触应力（MPa）；We为弹性势能密度函数，表示为刚度参数B，弹性体应变εve，弹性剪应变εse，温度T和胶结应力Cs的函数；Ke为材料的体积模量；$ \varepsilon _{ij}^p $为塑性应变；TGG为颗粒温度；$ e_{ij}^e $为弹性偏应变；ρs为固体颗粒相密度(kg/m3)；Ɍ为胶结材料影响系数，文中指水合物胶结材料；φ为胶结材料含量；C1，C2，C3，C4，C5为模型计算参数。
在非饱和土热力学模型基础上[21-22]，通过将理论模型中刚度参数B和胶结应力参数Cs引入温-压耦合系数δ(Δp, ΔT)，并在剪胀方程中考虑水合物饱和度的影响，最终建立了温度和孔隙压力耦合影响下水合物沉积物的热力学模型。
1.1 温-压耦合系数
水合物沉积物的试验结果表明温度越低，孔隙压力越大（试验反压），人工制备的水合物沉积物试样强度越高[10-11]。沉积物中颗粒与颗粒之间的胶结是通过甲烷气体在低温和高孔隙压力条件下结合孔隙水形成的[13]。在制备人工合成的水合物沉积物试样时，孔隙中仍存在未合成水合物的残余液态水和甲烷气体。该部分液体和气体在低温和高孔隙压力下一部分继续生成水合物，另一部分以冰的形式存在，最终使得沉积物的强度和刚度得到提高[16]。
将Teymouri等[24]提出的水合物沉积物状态参数δ(Δp, ΔT)引入水合物沉积物的热动力学理论模型中。水合物沉积物温-压耦合系数δ(Δp, ΔT)可以表示为当前状态与相平衡状态的相对距离，如图 1所示。该参数计算方式简单，物理含义明确，且能够考虑水合物沉积物形成过程中温度和孔隙压力的综合影响：
$$ \delta=\sqrt{\left[\delta_{\mathrm{T}}\left(T-T_{\mathrm{e}}\right)\right]^2+\left[\delta_{\mathrm{P}}\left(P-P_{\mathrm{e}}\right)\right]^2} 。 $$ (4) ![]()
图 1 水合物沉积物温-压耦合系数Figure 1. Temperature-pressure coupling coefficient of hydrate- bearing sediments式中：δT，δP分别为单位参数，δT=1 K-1，δP=1 MPa-1；T，P分别为当前状态下的温度（K）和孔隙压力（MPa）；Te，Pe分别为利用相平衡曲线计算出当前状态下对应的相平衡温度（K）和孔隙压力（MPa）。
相平衡线是描述处于相平衡状态的孔隙压力Pe和温度Te关系的曲线，如图 1所示。参考Gupta等[25]提出的相平衡公式，其具体表达式为
$$ {P_e} = {e_1}\exp \left( {{e_2} - \frac{{{e_3}}}{{{T_{\text{e}}}}}} \right) 。 $$ (5) 式中：e1，e2，e3分别为相平衡曲线参数，其中e1=1000 MPa，e2=38.98，e3=8533.8 K。
从图 1中可以看出，温-压耦合系数δ的平方值为温度和孔隙压力坐标内，状态点至温度相变点和孔隙压力相变点距离的平方和。当水合物沉积物状态点越远离温-压相变线时（A2点），参数δ的数值越大，对应水合物沉积物状态越稳定，其强度和刚度越大。
1.2 初始胶结应力的计算
关于胶结应力的计算，Uchida等[19]基于Masui等[8]试验数据给出了经验性表达式。该表达式以水合物含量为胶结应力的唯一控制因素。同时，在胶结应力衰减方面，考虑了剪应变的影响。然而，随着研究的深入，文献[11]发现温度和孔隙压力对含水合物沉积物强度和刚度有重要影响。而在胶结应力衰减方面, 体应变的影响不能被忽略。综合以上分析，提出一个新的胶结应力表达式及其衰减模式。为了简化计算，将温-压耦合系数δ对胶结应力的影响采用线性方式表达：
$$ C_{\mathrm{s}}=\alpha\left(S_{\mathrm{h}}\right)^\beta[1+\theta \delta(\Delta P, \Delta T)] \varPsi \quad\left(T \leqslant T_{\mathrm{e}}, P_{\mathrm{e}} \leqslant P\right) , $$ (6) $$ \varPsi=\exp \left(-\lambda_{\mathrm{v}} \varepsilon_{\mathrm{v}}^{\mathrm{p}}-\lambda_{\mathrm{s}} \varepsilon_{\mathrm{s}}^{\mathrm{p}}\right) 。 $$ (7) 式中：Sh为水合物饱和度（%）；Cs为胶结应力，其中初始胶结应力Cs0可以表示为Cs(Sh0)；α，β为胶结应力参数；θ为温-压耦合系数影响参数；Ψ为胶结应力衰减参数；$ {\lambda _\text{v}} $，$ {\lambda _\text{s}} $分别为塑性体应变和剪应变衰减系数。
在标定初始胶结应力时，Yan等[26]提出通过确定不同水合物饱和度条件下沉积物的屈服面来计算初始胶结应力Cs0。这种方法仅适用于强胶结情况。Shen等[27]在对试验结果的研究中发现含水合物沉积物的峰值剪切强度与其初始胶结强度存在对应关系。在峰值剪切强度之前，假定水合物沉积物试样的变形由颗粒之间的挤压造成，沉积物变形以弹性变形为主。此时固态水合物并没有发生大面积的破碎。而在峰值强度之后，伴随着水合物以及水合物与砂颗粒表面胶结物质破碎量增大，沉积物发生较大的塑性变形。沉积物强度降低，出现软化现象。根据此假定，本文给出了含水合物沉积物初始胶结应力的标定方式。
1.3 弹性势能密度函数构建
在颗粒热动力学中，颗粒之间的接触应力$ {\pi _{ij}} $可以通过弹性势能密度函数We求得。Bai等[22]证明通过适当简化，颗粒之间的接触应力$ {\pi _{ij}} $等于沉积物试样内部的有效应力$ {\sigma _{ij}} $。这样建立了弹性势能密度函数同有效应力之间的联系（式（1））。在构建非饱和土体的弹性势能密度函数中，Bai等[22]采用弹性应变和温度梯度的函数形式。该函数形式在验证非饱和土体的力学性质和饱和土体的热固结特性方面被证明是有效的。基于非饱和土体的弹性势能密度函数形式，通过引入胶结应力参数建立了水合物沉积物的弹性势能密度函数：
$$ W_{\mathrm{e}}=B\left[\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}\right)^{2.5}+\xi\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}\right)^{0.5}\left(\varepsilon_{\mathrm{s}}^{\mathrm{e}}\right)^2\right]+C_{\mathrm{s}} \varepsilon_{\mathrm{v}}^{\mathrm{e}} 。 $$ (8) 式中：ξ为模型中与材料摩擦性相关系数。
Bai等[22]提出参数B表征的是材料的胶结程度，B越高表示材料本征胶结强度越大。同时，参数B的衰减与孔隙比e的变化相关。在水合物沉积物中, 由于较低初始温度和较高孔隙压力的影响使得沉积物土体的强度得到增强，因此有必要建立初始刚度参数B0同温-压耦合系数δ关系。在非饱和土体刚度参数B的构建基础上提出将水合物沉积物的刚度参数表示为
$$ B = ({B_0} + {B_\delta })\exp \left( { - \frac{e}{\lambda }} \right) , $$ (9) $$ {B_\delta } = {B_0}\gamma {\text{ln(}}\delta + 1{\text{) (}}\delta \geqslant 0{\text{)}} 。 $$ (10) 式中：B0为水合物饱和度Sh=0 %条件下的刚度数值；Bδ表示温-压耦合系数对B0的影响；e为材料总孔隙比(e = ρLGs/ρ-1)；λ为材料等向压缩曲线斜率；Gs为材料的颗粒相对质量密度；ρ为混合材料物密度(kg/m3)；γ为温-压耦合系数对刚度参数的影响系数。
此外, 水合物沉积物由于内部胶结产生的结构性变化, 同重塑土在孔隙结构上存在差异。文献[28]中关于结构性土孔隙率的定义，引入有效孔隙率eeff概念, 建立水合物沉积物的孔隙率表达式：
$$ e=e_{\mathrm{eff}}+e_{\mathrm{h}}。 $$ (11) 式中：eeff为有效孔隙率(eeff =Vv/Vs)，表征的是在不考虑结构性条件下沉积物中的孔隙率；eh为水合物孔隙率(eh =She)，表征的是沉积物中由于胶结物质存在占据的空隙体积。
进一步化简后，水合物饱和度同沉积物土体内部总孔隙率e和有效孔隙率eeff的关系采用微分形式表达为
$$ d_{\mathrm{t}} e_{\mathrm{eff}}=\left(1-S_{\mathrm{h}}\right) d_{\mathrm{t}} e-\frac{e_{\mathrm{eff}}}{1-S_{\mathrm{h}}} d_{\mathrm{t}} S_{\mathrm{h}} 。 $$ (12) 将式（11），（12）代入式（9）中, 经过化简可以建立水合物沉积物刚度参数、弹性体积模量Ke和弹性剪切模量Ge表达式[21]：
$$ B = {B_0}\left[ {1 + \gamma {\text{ln}}(\delta + 1)} \right]\exp \left[ { - \frac{{{e_\text{eff}}}}{{\lambda (1 - {S_{\text{h}}})}}} \right] , $$ (13) $$ G_{\mathrm{e}}=B_0 \xi[1+\gamma \ln (\delta+1)]\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}\right)^{0.5} \exp \left[-\frac{e_{\mathrm{eff}}}{\lambda\left(1-S_{\mathrm{h}}\right)}\right] , $$ (14) $$ K_{\mathrm{e}}=B_0[1+\gamma \ln (\delta+1)] \exp \left[-\frac{e_{\mathrm{eff}}}{\lambda\left(1-S_{\mathrm{h}}\right)}\right] \cdot $$ $$ \left[2.5\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}\right)^{0.5}+0.5 \xi\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}\right)^{-1.5}\left(\varepsilon_{\mathrm{s}}^{\mathrm{e}}\right)^2\right] 。 $$ (15) 结合考虑式（1），（13）~（15）和土颗粒相的质量守恒方程[22]，建立含水合物沉积物的平均应力p′和偏应力q表达式：
$$ \begin{gathered} p^{\prime}=K_{\mathrm{e}} \varepsilon_{\mathrm{v}}^{\mathrm{e}}+C_{\mathrm{s}}=B_0[1+\gamma \ln (\delta+1)] \exp \left[-\frac{e_{\mathrm{eff}}}{\lambda\left(1-S_{\mathrm{h}}\right)}\right] . \\ {\left[2.5\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}\right)^{1.5}+0.5 \xi\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}\right)^{-0.5}\left(\varepsilon_{\mathrm{s}}^{\mathrm{e}}\right)^2\right]+C_{\mathrm{s}}, } \end{gathered} $$ (16) $$ q=\sqrt{\frac{3}{2} \pi_{k k} \pi_{k k}}=\sqrt{\frac{3}{2}\left(2 G_{\mathrm{e}} e_{i j}^{\mathrm{e}}\right)\left(2 G_{\mathrm{e}} e_{i j}^{\mathrm{e}}\right)} \\ =\sqrt{6} B_0 \xi[1+\gamma \ln (\delta+1)] \exp \left[-\frac{e_{\text {eff }}}{\lambda\left(1-S_{\mathrm{h}}\right)}\right]\left(\varepsilon_{\mathrm{v}}^{\mathrm{e}}\right)^{0.5} \varepsilon_{\mathrm{s}}^{\mathrm{e}} 。 $$ (17) 1.4 水合物沉积物剪胀方程
在构建水合物沉积物剪胀方程时，通过对Hyodo试验点的分析[27, 29]，假定不同水合物饱和度条件下沉积物的临界状态线斜率λc保持一致。同时，在临界状态条件下忽略温-压耦合系数的影响，只考虑水合物饱和度的影响。不同水合物饱和度条件下临界状态线参数eξ(Sh)不同，如图 2所示。将临界状态线参数eξ(Sh)设为[27]
$$ e_{\xi}\left(S_{\mathrm{h}}\right)=e_{\xi}+a\left(S_{\mathrm{h}}\right)^b 。 $$ (18) ![]()
图 2 水合物饱和度与临界孔隙比参数关系Figure 2. Relationship between hydrate saturation and critical void ratio parameters式中：eξ为水合物饱和度Sh=0条件下的临界状态线参数；а，b为临界状态线系数。
Been等[30]定义状态变量ψ用于描述颗粒物质密实程度。通过引入式（18）得到水合物饱和度影响下沉积物状态变量ψ(Sh)表达式:
$$ \psi\left(S_{\mathrm{h}}\right)=e-e_{\xi}-a\left(S_{\mathrm{h}}\right)^{\mathrm{b}}+\lambda_{\mathrm{c}} \ln \left(\frac{p^{\prime}}{p_{\mathrm{a}}}\right)。 $$ (19) 式中：$ {\lambda _\text{c}} $为水合物沉积物临界状态线参数，假定与水合物饱和度含量无关；pa为大气压力。
另一方面，水合物饱和度与沉积物临界应力比M(Sh)关系式可表示为[27]
$$ M({S_{\text{h}}}) = M + A{{\text{(}}{S_{\text{h}}})^B} 。 $$ (20) 式中：M为水合物饱和度Sh=0条件下的临界应力比；A，B为临界应力比系数。
参考Li等[31]建立的砂土剪胀方程, 将含水合物沉积物发生“相变状态”时的应力比Md定义为Md[M(Sh)，ψ(Sh), m]。引入式（19），（20），建立了水合物饱和度影响下沉积物的剪胀方程:
$$ d=d_0\left[\mathrm{e}^{m \mu\left(S_{\mathrm{h}}\right)}-\frac{\eta^*}{M\left(S_{\mathrm{h}}\right)}\right] 。 $$ (21) 式中：d0为水合物沉积物剪胀方程参数；η*为当前状态下的应力比；m为水合物沉积物剪胀方程参数。
通过剪胀比的定义($ d=d_{\mathrm{t}} \varepsilon_{\mathrm{v}}^{\mathrm{p}} / d_{\mathrm{t}} \varepsilon_{\mathrm{s}}^{\mathrm{p}}$)并结合考虑式（2），（21），得到水合物沉积物热力学本构关系中参数C1的表达式为
$$ C_1=d_0\left[\mathrm{e}^{m \mu /\left(S_{\mathrm{h}}\right)}-\frac{\eta^*}{M\left(S_{\mathrm{h}}\right)}\right]\left[\frac{\varepsilon_{\mathrm{s}}^{\mathrm{e}}}{\varepsilon_{\mathrm{v}}^{\mathrm{e}}+\frac{C_{\mathrm{s}}}{K_{\mathrm{e}}}}\right] 。 $$ (22) 2. 参数标定方法
计算参数共分为颗粒重组理论参数, 胶结应力参数和剪胀方程参数3类。其中颗粒重组理论参数在文献[22]中得到具体阐述，本文不再详细说明。值得注意的是颗粒重组理论中刚度参数通过回弹曲线斜率κ计算。由于试验数据的缺乏，在数值模拟验证时忽略了温-压耦合系数的影响（式（10）：$ \gamma $=0）。在参数敏感性分析时讨论了不同参数$ \gamma $数值下的沉积物力学性质变化。
对于胶结应力参数计算，基于Hyodo等[29]试验数据，计算各数据点对应的温-压耦合系数。将不同围压, 温-压耦合系数δ和水合物饱和度的数据点在p′-qmax平面内绘出。通过拟合方式得到对应水合物饱和度和温-压耦合系数的直线，如图 3所示。值得注意的是温-压耦合系数改变拟合直线的斜率，水合物饱和度的变化体现在截距的不同。在p′-qmax平面内，qmax=0时的p′即为对应沉积物试样的初始胶结应力。将得到的初始胶结应力数据在Sh-δ-Cs0三维空间内汇出, 通过拟合式（6）确定对应系数数值，如图 4所示。
![]()
图 3 沉积物峰值强度与水合物饱和度和温-压耦合系数关系Figure 3. Relationship among peak shear strength, hydrate saturation and temperature-pressure coupling coefficient![]()
图 4 初始胶结强度与水合物饱和度和温-压耦合系数关系Figure 4. Relationship among initial cementation strength, hydrate saturation and temperature-pressure coupling coefficient剪胀方程参数共分为两类，分别是水合物饱和度Sh=0时的参数M，eξ，d0和m和考虑水合物饱和度影响下的参数a，b，A，B。前者的计算方法同砂土剪胀参数计算方法一致[31]，文中不再详细叙述。对于后者，系数(a, b)和(A, B)分别利用不同水合物饱和度条件下的临界空隙比eξ(Sh)和临界应力比M(Sh)数值，通过式（18），（20）标定得到[27]，如图 2，5所示。
![]()
图 5 水合物饱和度与临界应力比关系Figure 5. Relationship between hydrate saturation and critical stress ratio3. 模型验证
以Hyodo等[32]试验数据为参考点，通过对本构方程（式（1）~（3），（15），（16），（22））利用Matlab 2017进行编程，采用龙格-库塔迭代方法，将数值计算结果与试验结果进行对比。其中，Hyodo试验以Toyoura砂为土骨架材料（Dr=90%），利用气饱和制样方式对不同孔隙压力（5，10，12 MPa）、温度（1℃，5℃）、围压（1，3 MPa）和水合物饱和度（Sh= 23.2%~47.7%）的试样进行三轴排水剪切试验。对应的颗粒热动力学理论参数：B0 =6.03×107 MPa，$ \lambda $=0.13，ξ=6，C3/s2=10×104，C5=9×102 kg/m3/s，$ \gamma $=0，Gs=2.65。颗粒热动力学理论参数：$ {\lambda _{\text{v}}} $=10，$ {\lambda _{\text{s}}} $=6，$ \alpha $=0.59，β=0.73，θ=0.04。剪胀方程参数：M=1.16，eξ=1.13，d0=1.6，m=1.3，a=0.50，b=1.78，A=0.36，B=1.20，$ {\lambda _{\text{c}}} $=1.126。
3.1 围压对沉积物力学性质影响
三轴排水剪切试验中相同水合物饱和度和温-压耦合系数，不同围压条件下各点的应力-应变与体应变关系如图 6所示。其中，空心点表示试验数据，实线是模拟数据。从模拟结果可以看出，基于颗粒重组理论的水合物沉积物模型能够很好的体现围压对沉积物力学性质的影响，主要体现在两个方面：①在应力-应变关系中，峰值偏应力和初始阶段曲线斜率随着围压的增大而增大，说明围压与沉积物试样的强度与弹性模量成正相关。此外，在剪切阶段的后期（15%应变）模拟和试验结果均呈现出软化现象。②在应变-体应变关系中，剪胀效应随着围压的增大逐渐减弱, 这与文献[28，29]中得出的结论一致。
![]()
图 6 不同围压下模拟值与试验值对比Figure 6. Comparative results by proposed model and tests under different confining pressures3.2 水合物饱和度对沉积物力学性质影响
图 7（a），（b）描述了在相同围压和温-压耦合系数条件下，不同水合物饱和度沉积物试样轴向应变-偏应力/体应变的模型预测结果。从模拟结果与试验结果对比分析可以看出，模型能够很好的描述水合物饱和度对沉积物试样的强度、初始弹性模型和剪胀性的增强作用。同时，随着水合物饱和度的增大，沉积物试样逐步表现出应变软化性质。
![]()
图 7 不同水合物饱和度下模拟值与试验值对比Figure 7. Comparative results by proposed model and tests under different hydrate saturations3.3 温-压耦合系数对沉积物力学性质影响
图 8考虑了相同孔隙压力（10 MPa），不同剪切温度（T为1℃，5℃）的影响，其对应的温-压耦合系数δ为14.54，10.39。从图 8结果可以看出建立的模型能够反映不同温-压耦合系数下沉积物试样的应力和应变演化特性，具体表现出随着参数δ数值的增大，沉积物试样的强度、初始弹性模量和剪胀特性增大。该现象可以理解为表征着赋存环境影响的温-压耦合系数增大提高了微观层次土颗粒之间的胶结作用。一方面使得水合物沉积物表现出更大的峰值强度和初始弹性模量。另一方面，从颗粒重组理论角度考虑，破碎后的水合物（应变软化阶段）填充进土骨架的孔隙中，进一步提高了颗粒重组性能，在宏观上表现为剪胀性能提高。
![]()
图 8 不同温-压耦合系数(温度引起)下模拟值与试验值对比Figure 8. Comparative results by proposed model and tests under different temperature-pressure coupling coefficients caused by the temperatures但在局部区域，如三轴排水剪切的初始阶段，理论模型低估了水合物沉积物试样的强度。造成差异的原因是本文在对刚度参数B计算中忽略了温-压耦合系数的影响，从而低估了模型中刚度参数B0的数值。这一点在参数敏感性分析时进行详细介绍。
图 9呈现了相同剪切温度（T=5℃），不同孔隙压力（5，10 MPa）时（δ为2.21，10.39），水合物沉积物试样的轴向应变-偏应力/体应变的模拟计算结果图。通过与试验结果对比发现模型能够预测沉积物试样因温-压耦合系数变化产生的物理力学性质不同。由于低估了刚度参数B0的数值，在三轴排水剪切的初始阶段，模型结果的预测值与实验值存在偏差。
![]()
图 9 不同温-压耦合系数(孔隙压力引起)下模拟值与试验值对比Figure 9. Comparative results by proposed model and tests under different temperature-pressure coupling coefficients caused by the water pressures4. 参数敏感性分析
从图 8，9的模拟结果和试验结果对比可以看出由于计算初始刚度参数B0数值时忽略了温-压耦合系数的影响，使得B0的数值被低估，最终在三轴排水剪切初始阶段造成模拟结果与试验结果差异。从式（10）中可以看出，温-压耦合系数对参数B0的影响通过系数$ \gamma $体现。在图 10中分别考察了4种$ \gamma $数值条件下（$ \gamma $为0，0.2，0.5，1.0），沉积物试样的力学性质变化。其中水合物饱和度Sh=26.6%，围压3 MPa，孔隙压力10 MPa。
![]()
图 10 刚度系数γ对沉积物力学性质影响Figure 10. Effects of stiffness coefficient on mechanical properties of sedimentsBai等[22]将参数B0定义为土颗粒材料的胶结程度。参数B0越大，表明材料自身的强度和刚度越大，这一点同超固结效应类似。因此，在外力作用下试样的初始弹性模量和峰值强度随着系数$ \gamma $的增大而增强（图 10（a））。值得注意的是不同系数$ \gamma $条件下试样的破坏强度（轴向应变为15%）接近。该现象说明了参数B0的增大是由于微观上颗粒之间胶结作用引起的。在达到峰值强度后，随着土颗粒之间胶结效应的衰减（式（7）），土骨架的影响逐渐成为主导因素。另一方面，在应变软化阶段随着水合物胶结作用的衰减，破碎后的水合物填充进土骨架孔隙中。从颗粒重组理论角度分析，破碎后的水合物的填充作用进一步增强了颗粒之间的重组效应，宏观表现为剪胀性提高。
另一方面，本文提出了两种因素影响下水合物沉积物胶结应力衰减因子，即体应变衰减因子$ {\lambda _{\text{v}}} $和剪应变衰减因子$ {\lambda _{\text{s}}} $。文献[26]对剪应变衰减因子$ {\lambda _{\text{s}}} $的参数影响性分析进行了集中的讨论，文中不再详述。图 11则是考虑3种体应变衰减因子（$ {\lambda _{\text{v}}} $为1，10，20，30）对沉积物试样力学性质的影响。在图 11（a）中，$ {\lambda _{\text{v}}} $增大，沉积物试样的峰值强度和初始弹性模量没有明显变化，而应变软化性质增强。$ {\lambda _{\text{v}}} $表征的是三轴排水剪切过程中沉积物内部胶结应力的衰减程度，$ {\lambda _{\text{v}}} $越大，表明胶结应力衰减越迅速，宏观表现是应变软化性质明显。模拟结果与预期值相符。此外，随着$ {\lambda _{\text{v}}} $的增大，沉积物试样的剪胀特性明显增强。基于颗粒热动力学理论分析，随着胶结应力的快速衰退，沉积物内部以胶结形式连接的土骨架坍塌，最终使得内部土颗粒产生重组。而$ {\lambda _{\text{v}}} $越大，衰退效应越明显，土颗粒重组效应越强，宏观表现是材料的剪胀性能增强。
![]()
图 11 体应变衰减系数λv对沉积物力学性质影响Figure 11. Effects of degradation coefficient of volumetric strain on mechanical properties of sediments5. 结论
基于颗粒重组理论，建立了含水合物沉积物的温度-孔隙压力-应力耦合热力学模型。通过对含水合物沉积物进行三轴排水剪切试验模拟，并将模拟结果与Hyodo的室内试验结果对比，验证了模型的有效性。此外，分析了沉积物剪切时的温度和孔隙压力变化对力学性质的影响规律，并对模型建立过程中的刚度参数$ \gamma $和胶结应力衰减参数$ {\lambda _{\text{v}}} $进行敏感性分析。
（1）引入温-压耦合系数和考虑水合物饱和度影响下剪胀方程，建立含水合物沉积物热力学模型能够有效的描述沉积物试样的应力-应变和体应变的变化规律同温-压耦合系数，围压和水合物饱和度的变化关系。
（2）从微观角度分析，降低水合物沉积物赋存时的温度或者升高孔隙压力能够提高沉积物内部初始胶结强度。对应的宏观力学响应表现为峰值剪切强度、初始弹性模量、应变软化性能和剪胀性得到提高。
（3）模型中的刚度参数$ \gamma $与沉积物的峰值剪切强度和初始弹性模量呈正相关，但对沉积物的应变软化性能和残余强度影响较小。而衰减参数$ {\lambda _{\text{v}}} $越大，沉积物试样的应变软化性越明显，残余强度越低，但对峰值强度和弹性模量影响较弱。
（4）沉积物试样赋存时的温-压耦合系数对峰值剪切强度、弹性模量和应变软化性影响明显。而剪胀性受温-压耦合系数影响有限。该系数可以用来评价复杂温-压环境下水合物沉积物的力学特性。
-
![]()
图 1 水合物沉积物温-压耦合系数
Figure 1. Temperature-pressure coupling coefficient of hydrate- bearing sediments
![]()
图 2 水合物饱和度与临界孔隙比参数关系
Figure 2. Relationship between hydrate saturation and critical void ratio parameters
![]()
图 3 沉积物峰值强度与水合物饱和度和温-压耦合系数关系
Figure 3. Relationship among peak shear strength, hydrate saturation and temperature-pressure coupling coefficient
![]()
图 4 初始胶结强度与水合物饱和度和温-压耦合系数关系
Figure 4. Relationship among initial cementation strength, hydrate saturation and temperature-pressure coupling coefficient
![]()
图 5 水合物饱和度与临界应力比关系
Figure 5. Relationship between hydrate saturation and critical stress ratio
![]()
图 6 不同围压下模拟值与试验值对比
Figure 6. Comparative results by proposed model and tests under different confining pressures
![]()
图 7 不同水合物饱和度下模拟值与试验值对比
Figure 7. Comparative results by proposed model and tests under different hydrate saturations
![]()
图 8 不同温-压耦合系数(温度引起)下模拟值与试验值对比
Figure 8. Comparative results by proposed model and tests under different temperature-pressure coupling coefficients caused by the temperatures
![]()
图 9 不同温-压耦合系数(孔隙压力引起)下模拟值与试验值对比
Figure 9. Comparative results by proposed model and tests under different temperature-pressure coupling coefficients caused by the water pressures
![]()
图 10 刚度系数γ对沉积物力学性质影响
Figure 10. Effects of stiffness coefficient on mechanical properties of sediments
-
[1] INADA N, YAMAMOTO K. Data report: hybrid pressure coring system tool review and summary of recovery result from gas-hydrate related coring in the Nankai Project[J]. Marine and Petroleum Geology, 2015, 66: 323-345. doi: 10.1016/j.marpetgeo.2015.02.023
[2] MAKOGON Y F. Natural gas hydrates—a promising source of energy[J]. Journal of Natural Gas Science and Engineering, 2010, 2(1): 49-59. doi: 10.1016/j.jngse.2009.12.004
[3] 颜荣涛, 韦昌富, 魏厚振, 等. 水合物形成对含水合物砂土强度影响[J]. 岩土工程学报, 2012, 34(7): 1234-1240. http://cge.nhri.cn/cn/article/id/14631 YAN Rongtao, WEI Changfu, WEI Houzhen, et al. Effect of hydrate formation on mechanical strength of hydrate-bearing sand[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(7): 1234-1240. (in Chinese) http://cge.nhri.cn/cn/article/id/14631
[4] 刘昌岭, 李彦龙, 孙建业, 等. 天然气水合物试采: 从实验模拟到场地实施[J]. 海洋地质与第四纪地质, 2017, 37(5): 12-26. https://www.cnki.com.cn/Article/CJFDTOTAL-HYDZ201805001.htm LIU Changling, LI Yanlong, SUN Jianye, et al. Gas hydrate production test: from experimental simulation to field practice[J]. Marine Geology & Quaternary Geology, 2017, 37(5): 12-26. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-HYDZ201805001.htm
[5] 王维希, 张春生, 吴颜雄, 等. 联合深海地热开采天然气水合物技术展望[J]. 现代化工, 2021, 41(9): 17-21. https://www.cnki.com.cn/Article/CJFDTOTAL-XDHG202109004.htm WANG Weixi, ZHANG Chunsheng, WU Yanxiong, et al. Prospects on technology for combining deep-sea geothermal energy with exploitation of natural gas hydrate[J]. Modern Chemical Industry, 2021, 41(9): 17-21. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-XDHG202109004.htm
[6] WU P, LI Y H, LIU W G, SUN X, KONG X J, SONG Y C. Cementation failure behavior of consolidated gas hydrate-bearing sand[J]. Journal of Geophysical Research: Solid Earth, 2020, 125(1): e2019JB018623. doi: 10.1029/2019JB018623
[7] LUO T T, LI Y H, MADHUSUDHAN B N, et al. Deformation behaviors of hydrate-bearing silty sediment induced by depressurization and thermal recovery[J]. Applied Energy, 2020, 276: 115468. doi: 10.1016/j.apenergy.2020.115468
[8] MASUI A, MIYAZAKI K, HANEDA H, OGATA Y, AOKI K. Mechanical properties of natural gas hydrate bearing sediments retrieved from eastern Nankai trough[C]//Offshore Technology Conference. Houston, 2008.
[9] LI Y H, SONG Y C, LIU W G, et al. Analysis of mechanical properties and strength criteria of methane hydrate-bearing sediments[J]. International Journal of Offshore and Polar Engineering, 2012, 22: 290-296.
[10] SHEN S, SUN X, WANG L, et al. Effect of temperature on the mechanical properties of hydrate-bearing sand under different confining pressures[J]. Energy & Fuels, 2021, 35: 4106-4117.
[11] HYODO M, NAKATA Y, YOSHIMOTO N, et al. Basic research on the mechanical behavior of methane hydrate-sediments mixture[J]. Soils and Foundations, 2005, 45(1): 75-85.
[12] WINTERS W J, WAITE W F, MASON D H, et al. Methane gas hydrate effect on sediment acoustic and strength properties[J]. Journal of Petroleum Science and Engineering, 2007, 56(1/2/3): 127-135.
[13] CHAOUACHI M, FALENTY A, SELL K, et al. Microstructural evolution of gas hydrates in sedimentary matrices observed with synchrotron X-ray computed tomographic microscopy[J]. Geochemistry, Geophysics, Geosystems, 2015, 16(6): 1711-1722. doi: 10.1002/2015GC005811
[14] TOKUNAGA T. Physicochemical controls on adsorbed water film thickness in unsaturated geological media[J]. Water Resources Research, 2011, 47(8): W08514.
[15] 蒋明镜, 朱方园. 不同温压环境下深海能源土力学特性离散元分析[J]. 岩土工程学报, 2014, 36(10): 1761-1769. doi: 10.11779/CJGE201410001 JIANG Mingjing, ZHU Fangyuan. DEM investigation on mechanical properties of methane hydrate bearing soils under different temperatures and pore-water pressures[J]. Chinese Journal of Geotechnical Engineering, 2014, 36(10): 1761-1769. (in Chinese) doi: 10.11779/CJGE201410001
[16] 颜荣涛, 张炳晖, 杨德欢, 等. 不同温-压条件下含水合物沉积物的损伤本构关系[J]. 岩土力学, 2018, 39(12): 4421-4428. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201812016.htm YAN Rongtao, ZHANG Binghui, YANG Dehuan, et al. Damage constitutive model for hydrate-bearing sediment under different temperature and pore pressure conditions[J]. Rock and Soil Mechanics, 2018, 39(12): 4421-4428. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201812016.htm
[17] 刘乐乐, 张旭辉, 刘昌岭, 等. 含水合物沉积物三轴剪切试验与损伤统计分析[J]. 力学学报, 2016, 48(3): 720-729. https://www.cnki.com.cn/Article/CJFDTOTAL-LXXB201603021.htm LIU Lele, ZHANG Xuhui, LIU Changling, et al. Triaxial shear tests and statistical analyses of damage for methane hydrate-bearing sediments[J]. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(3): 720-729. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-LXXB201603021.htm
[18] YAMAGUCHI T, AOKI K, TENMA N, et al. A nonlinear elastic model for triaxial compressive properties of artificial methane-hydrate-bearing sediment samples[J]. Energies, 2012, 5(10): 4057-4075. doi: 10.3390/en5104057
[19] UCHIDA S, SOGA K, YAMAMOTO K. Critical state soil constitutive model for methane hydrate soil[J]. Journal of Geophysical Research: Solid Earth, 2012, 117: B03209.
[20] SÁNCHEZ M, GAI X R, SANTAMARINA J C. A constitutive mechanical model for gas hydrate bearing sediments incorporating inelastic mechanisms[J]. Computers and Geotechnics, 2017, 84: 28-46. doi: 10.1016/j.compgeo.2016.11.012
[21] YANG G C, BAI B, LIU Y, et al. Constitutive modeling for undrained shear behavior of gassy sand considering energy dissipation at the mesoscopic level[J]. Ocean Engineering, 2021, 219: 108307. doi: 10.1016/j.oceaneng.2020.108307
[22] BAI B, ZHOU R, CAI G Q, et al. Coupled thermo- hydro-mechanical mechanism in view of the soil particle rearrangement of granular thermodynamics[J]. Computers and Geotechnics, 2021, 137: 104272. doi: 10.1016/j.compgeo.2021.104272
[23] JIANG Y M, LIU M. Granular solid hydrodynamics[J]. Granular Matter, 2009, 11(3): 139-156. doi: 10.1007/s10035-009-0137-3
[24] TEYMOURI M, SÁNCHEZ M, SANTAMARINA J C. A pseudo-kinetic model to simulate phase changes in gas hydrate bearing sediments[J]. Marine and Petroleum Geology, 2020, 120: 104519. doi: 10.1016/j.marpetgeo.2020.104519
[25] GUPTA S, HELMIG R, WOHLMUTH B. Non-isothermal, multi-phase, multi-component flows through deformable methane hydrate reservoirs[J]. Computational Geosciences, 2015, 19(5): 1063-1088. doi: 10.1007/s10596-015-9520-9
[26] YAN R T, WEI C F. Constitutive model for gas hydrate-bearing soils considering hydrate occurrence habits[J]. International Journal of Geomechanics, 2017, 17(8): 04017032. doi: 10.1061/(ASCE)GM.1943-5622.0000914
[27] SHEN J, CHIU C F, NG C W W, et al. A state-dependent critical state model for methane hydrate-bearing sand[J]. Computers and Geotechnics, 2016, 75: 1-11. doi: 10.1016/j.compgeo.2016.01.013
[28] DE LA FUENTE M, VAUNAT J, MARÍN-MORENO H. A densification mechanism to model the mechanical effect of methane hydrates in sandy sediments[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2020, 44(6): 782-802. doi: 10.1002/nag.3038
[29] HYODO M, YONEDA J, YOSHIMOTO N, et al. Mechanical and dissociation properties of methane hydrate-bearing sand in deep seabed[J]. Soils and Foundations, 2013, 53(2): 299-314. doi: 10.1016/j.sandf.2013.02.010
[30] BEEN K, JEFFERIES M G. A state parameter for sands[J]. Geotechnique, 1985, 35(2): 99-112. doi: 10.1680/geot.1985.35.2.99
[31] LI X S, DAFALIAS Y F. Dilatancy for cohesionless soils[J]. Géotechnique, 2000, 50(4): 449-460. doi: 10.1680/geot.2000.50.4.449
[32] HYODO M, LI Y H, YONEDA J, et al. Mechanical behavior of gas-saturated methane hydrate-bearing sediments[J]. Journal of Geophysical Research: Solid Earth, 2013, 118: 5185-5194. doi: 10.1002/2013JB010233
-
其他相关附件
下载:
计量
- 文章访问数: 250
- HTML全文浏览量: 74
- PDF下载量: 50
