[发明专利]一种考虑水合物分解的沉积物多场耦合模型的建模方法

权利要求书

1.一种考虑水合物分解的沉积物多场耦合模型的建立方法，其特征在于，步骤如下：

(一)基于孔隙介质力学推导多场耦合控制方程

根据孔隙介质力学，α相物质体积分数为:

其中，α为s、w、g或h，分别代表沉积物颗粒、水、气体和水合物；V表示总体积，V_α表示α相物质的体积；水合物、孔隙水和孔隙气饱和度分别定义为s_h＝n_h/n、s_w＝n_w/n和s_g＝n_g/n，其中，n是孔隙率；三者之间关系如下所示：

s_w+s_h+s_g＝1 (2)

根据孔隙介质力学，α相物质的质量守恒方程表示为：

其中，为α相物质速度散度，ρ_α表示α相物质的密度，dm_α/dt表示α相物质的质量累积，即源或汇；

沉积物颗粒的连续性方程为：

以沉积物颗粒作为骨架，则有ε_v为体积应变；

对于水和气体，其连续性方程为：

水合物连续性方程为：

沉积物骨架由沉积物颗粒与水合物共同构成，且不考虑水合物从沉积物骨架中挤出，故v_s＝v_h；

方程中的质量累积项满足化学反应过程中的质量守恒，因此

将流速使用达西流速表示v_w＝q_w/n_w+v_s和v_g＝q_g/n_g+v_s分别代入到连续性方程(6)和(7)中，其中，v_w、v_g和v_s分别是孔隙水、孔隙气以及沉积物骨架的速度，并且忽略n_w、n_g和密度的空间分布对质量守恒的影响，得：

通过孔隙介质力学理论描述水合物沉积物的能量守恒；满足机械能守恒定律，不考虑内力做功对温度变化的微弱影响；方程如下所示：

其中，T为温度，R_h为水合物反应速率，导热系数K_T利用叠加原理获得：K_T＝n_sK_Ts+n_gK_Tg+n_hK_Th+n_wK_Tw；K_Ts、K_Tg、K_Th和K_Tw分别为固体、气体、水合物和水的导热系数；c_T是水合物沉积物的总热容，满足c_T＝ρ_sn_sc_s+ρ_gn_gc_g+ρ_hn_hc_h+ρ_wn_wc_w，其中，c_s、c_g、c_h、c_w分别为沉积物颗粒、气体、水合物和水的比热；ΔH为每mol水合物相变引起的热量变化，通过Kamath回归方程计算，满足ΔH＝c₁+c₂T，c₁和c₂为两个回归参数；

力平衡方程的具体描述为：

其中，p_w为孔隙水压力，p_g为孔隙气压力，σ为总应力，f_α为α相物质流体的渗透阻力；根据有效应力原理，截面上应力等于有效应力与孔隙压力之和，表示为：

σ＝σ'-p_pδ (16)

其中，σ'为有效应力，p_p为孔隙压力；δ为单位矩阵；对于水合物沉积物，其孔隙压力表达式如下所示：

(二)COMSOL Multiphysics PDE建模

控制方程包括方程(10)-(13)，其中偏微分方程通过有限元离散化为常微分方程进行求解，采用的单元类型为混阶Lagrange单元；选取Lagrange插值形函数N_u，N_pw，N_pg，N_T，N_qw，N_qg，N_qT代入到控制方程中，并通过变分原理得到常微分方程：

K_uu:u+K_upgp_g+K_upwp_w＝F_u (18)

其中，

式中：L₁为微分算子矩阵，D为弹性矩阵，ε为应变，ε^p为塑性应变，t为面力张量，f为体力张量，q_T为热流，β_h，β_w和β_g分别为水合物、水以及气体的热膨胀系数，K_h为水合物沉积物的固有渗透系数，k_w和k_g分别为水和气体的相对渗透系数，μ_w和μ_g分别是水和气体的粘滞性系数，B_w为水的体积膨胀系数，M_g、M_h以及M_w分别为气体、水合物和水的摩尔质量，R_h为水合物反应速率，δ为单位矩阵，n_法为边界面法向矢量，g_重为重力加速度；

联立上述方程，对指定的工况设置初值和边值条件，得到整个求解初值和边值问题的方程列式；

求解常微分方程组(18)-(21)，首先进行时域离散，将常微分方程离散为代数方程系统来进行求解；常微分方程系统写成：

0＝M(U,t) (23)

其中，L、N_F、M分别代表刚度、约束力与约束，U为自由度；在求解该微分方程之前先化简拉格朗日乘子Λ；当0＝M(U,t)是线性的或时间相关时，且当约束力雅克比为常数时，消去约束M；其他情况下，则保留约束M进行求解；

对于非线性问题，采用牛顿迭代法，其线性化公式表示为：

NΔU＝M (25)

其中，为刚度矩阵，称为阻尼矩阵，为质量矩阵，ΔU为每一步计算得到的U的增量；

假设自由度U通过如下关系化为一组无量纲化后的自由度表示：

S是无量纲化系数矩阵；对线性方程组无量纲化后得到：

其中，并且式中E_r是一个对角矩阵以使为单位阵；

假设U_i是某一时间步上自由度i的解，E_i是求解器在该时间步上对U_i的绝对误差A_i做的估计，则表示为：

其中，E_rror为相对误差，M_f是参与耦合问题的物理场的个数，N_j是第j个物理场的自由度个数，该收敛准则等价于E_iA_i+E_r|U_i|；

(三)基于热力学方法的水合物沉积物临界状态本构模型的返回映射算法构造

从弹性力学中应力σ与应变ε关系出发，计算试应力σ_tri和修正应力σ_corr：

其中，C为热膨胀系数矢量，E为弹性模量；试应力为上一步传递过来的应力σ_n加上通过线弹性计算得到的应力增量：

修正应力是总的应力减去试应力得到的部分：

其中，h_硬是硬化参数矢量，下角标n+1表示第n+1步，C为热膨胀系数矢量，T为温度，E为弹性模量；弹性模量的变化由水合物对沉积物模量贡献给出：

其中，χ是结构因子，sh为水合物饱和度；

建立误差函数对误差函数进行离散：

f(σ,h_硬,sh_n+1)＝0 (35)

其中，a,b分别为塑性应变矢量误差，f为屈服面函数，sh_n+1为第n+1步水合物饱和度，λ为拉格朗日乘子，因为积分格式为隐式，所以第n+1步的水合物饱和度直接由质量守恒方程获得；r为塑性应变增量的方向，v为硬化参数矢量的方向；

通过离散得到线性代数方程，获得第k步的线性代数方程：

a^(k)-Δε^p(k)+δλ^(k)r^(k)+Δλ^(k)Δr^(k)＝0 (36)

b^(k)-Δh_硬^(k)+δλ^(k)v^(k)+Δλ^(k)Δv^(k)＝0 (37)

公式(31)代入到(36)-(38)中，得到：

b^(k)-Δh_硬^(k)+δλ^(k)v^(k)+Δλ^(k)Δv^(k)＝0 (40)

对上述方程进行重新整理，得到：

其中：

将应力和硬化参数增量代入f的表达式中解出δλ：

模型中，屈服面f的表达式为：

该模型引入应力比参数γ来控制干面和湿面的比例，由于受应力比的影响，屈服面将不再保持为椭圆，通过插值的方法定义转移应力ρ'；其中，M_临界为临界状态线斜率，α是另一个与应力间距比相关的参数；

p’为有效围压，q为广义剪应力，硬化参数矢量h_硬包括R、p′_cd、p′_cs以及p′_cc，分别描述屈服前的非线性弹性、水合物胶接结构变化引起的剪胀行为、沉积物的固结过程以及胶接结构破坏引起的黏聚力降低；

剪胀规律满足：

(四)将本构模型嵌入到COMSOL Multiphysics中实现水合物分解过程沉积物变形预测。