[發明專利]一種考慮水合物分解的沉積物多場耦合模型的建模方法

有效

• 申請號: 201710416098.2
• 申請日: 2017-06-06
• 公開（公告）號: CN107122571B
• 公開（公告）日: 2020-09-11
• 發明（設計）人: 宋永臣;李洋輝;孫翔;劉衛國;楊明軍;趙佳飛;劉瑜;張毅;王大勇;趙越超;蔣蘭蘭;凌錚;羅昊
• 申請（專利權）人: 大連理工大學
• 主分類號: G06F30/23
• 分類號: G06F30/23
• 代理公司: 大連理工大學專利中心 21200
• 代理人: 溫福雪;侯明遠

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中實現水合物分解過程沉積物變形預測。