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

摘要

本发明属于岩土工程技术领域，提供了一种考虑水合物分解的沉积物多场耦合模型的建模方法。本发明基于孔隙介质力学基本理论，结合水合物沉积物的力学特性、流动特性以及传热特性，给出了描述水合物沉积物物理力学行为的控制方程和本构方程。利用COMSOL Multiphysics软件中的PDE模块实现温度、压力以及饱和度的控制方程离散，采用结构力学模块中的外部材料接口，嵌入了基于热力学方法的水合物沉积物临界状态本构模型。同时，采用向后欧拉积分形式的返回映射算法，利用C++语言编写模型程序，实现考虑水合物分解过程的沉积物力学行为数值模拟与分析。

主权项

一种考虑水合物分解的沉积物多场耦合模型的建立方法，其特征在于，步骤如下：(一)基于孔隙介质力学推导多场耦合控制方程根据孔隙介质力学，α相物质体积分数为其中，α为s、w、g或h，分别代表沉积物颗粒、水、气体和水合物；V表示体积，Vα表示α相物质的体积；水合物、孔隙水和孔隙气饱和度分别定义为sh＝nh/n、sw＝nw/n和sg＝ng/n，其中，n是孔隙率；三者之间关系如下所示：sw+sh+sg＝1 (2)dswdt+dsgdt+dshdt=0---(3)]]>根据孔隙介质力学，α相物质的质量守恒方程表示为：d(ραnα)dt+ραnα▿·v=dmαdt---(4)]]>其中，为速度散度，ρα表示α相物质的密度，dmα/dt表示α相物质的质量累积，即源或汇；沉积物颗粒的连续性方程为：-11-ndndt+▿·vs=0---(5)]]>以沉积物颗粒作为骨架，则有εv为体积应变；对于水和气体，其连续性方程为：ndswdt+swdndt+nwρwdρwdt+nw▿·vw=1ρwdmwdt---(6)]]>ndsgdt+sgdndt+ngρgdρgdt+ng▿·vg=1ρgdmgdt---(7)]]>水合物连续性方程为：ndshdt+shdndt+nhρhdρhdt+nh▿·vh=1ρhdmhdt---(8)]]>沉积物骨架由沉积物颗粒与水合物共同构成，且不考虑水合物从沉积物骨架中挤出，故vs＝vh；方程中的质量累积项满足化学反应过程中的质量守恒，因此dmwdt+dmgdt+dmhdt=0---(9)]]>将流速使用达西流速表示vw＝qw/nw+vs和vg＝qg/ng+vs分别代入到连续性方程(6)和(7)中，其中，vw、vg和vs分别是孔隙水、孔隙气以及土骨架的速度，并且忽略nw、ng和密度的空间分布对质量守恒的影响，得：ndswdt+swdϵvdt+nwρwdρwdt+▿·qw=1ρwdmwdt---(10)]]>ndsgdt+sgdϵvdt+ngρgdρgdt+▿·qg=1ρgdmgdt---(11)]]>通过孔隙介质力学理论描述水合物沉积物的能量守恒；满足机械能守恒定律，不考虑内力做功对温度变化的微弱影响；方程如下所示：cTdTdt+▿·(cwρwqwT)+▿·(cgρgqgT)+▿·(-KT▿T)=-ΔHRh---(12)]]>其中，导热系数KT利用叠加原理获得：KT＝nsKTs+ngKTg+nhKTh+nwKTw；KTs、KTg、KTh和KTw分别为固体、气体、水合物和水的导热系数；cT是水合物沉积物的总热容，满足cT＝ρsnscs+ρgngcg+ρhnhch+ρwnwcw，其中，cs、cg、ch、cw分别为沉积物颗粒、气体、水合物和水的比热；ΔH为每mol水合物相变引起的热量变化，通过Kamath回归方程计算，满足ΔH＝c1+c2T，c1和c2为两个回归参数；力平衡方程的具体描述为：▿·(-σ)=(nsρs+nwρw+ngρg+nhρh)g---(13)]]>▿pw+fw=ρwg---(14)]]>▿pg+fg=ρgg---(15)]]>其中，σ为总应力，fα为α相物质流体的渗透阻力；根据有效应力原理，截面上应力等于有效应力与孔隙压力之和，表示为：σ＝σ'‑ppδ(16)其中，σ'为有效应力，pp为孔隙压力；δ为单位矩阵；对于水合物沉积物，其孔隙压力表达式如下所示：其中，pw和pg分别为孔隙水压力和孔隙气压力；(二)COMSOL Multiphysics PDE建模控制方程包括方程(8)、(10)‑(13)，其中偏微分方程通过有限元离散化为常微分方程进行求解，采用的单元类型为混阶Lagrange单元；选取Lagrange插值形函数Nu，Npw，Npg，NT，Nqw，Nqg，NqT代入到控制方程中，并通过变分原理得到常微分方程：Kuu:u+Kupgpg+Kupwpw＝Fu  (18)Cpwu·dudt+Cpwpwdpwdt+Kpwpwpw+Cpwpgdpgdt+KpwTdTdt=Fpw---(19)]]>Cpgu·dudt+Cpgpwdpwdt+Cpgpgdpgdt+Kpgpgpg+CpgTdTdt=Fpg---(20)]]>CTTdTdt+KTTT+KTTpwT+KTTpgT=FT---(21)]]>其中，Kupg=∫Ω▿·Nu·(-sgNpgsw+sgδ)dV;]]>Kupw=∫Ω▿·Nu·(-swNpwsw+sgδ)dV;]]>Fu=∫∂ΩNu·t·NtdS+∫ΩNu·f·NfdV;]]>Cpwu=∫ΩNpw(sw-sh∂sw∂sh)δ:L:NudV;]]>Cpwpw=∫ΩNpw((nwBw-n∂sw∂sedsedpc)Npw)dV;]]>Kpwpw=∫Ω▿Npw·Khkwμw(▿Npw)dV;]]>Cpwpg=∫ΩNpw(n∂sw∂sedsedpcNpg)dV;]]>KpwT=∫ΩNpw(-(nwβw-nhβh∂sw∂sh)NT)dV;]]>Fpw=∫Ω▿Npw·KhkwρwμwgdV+∫ΩNpw(61ρwMw+1ρh∂sw∂shMh)RhdV+Σi∫∂ΩNpwn·Nqw·qwdS;]]>Cpgu=∫ΩNpg(sg-sh∂sg∂sh)δ:L:NudV;]]>Cpgpw=∫Ω(-Npgn∂sg∂sedsedpcNpw)dV;]]>Cpgpg=∫ΩNpg((ngBg+n∂sg∂sedsedpc)Npg)dV;]]>Kpgpg=∫Ω▿Npg·Khkgμg(▿Npg)dV;]]>CpgT=∫ΩNpg(nhβh∂sg∂sh-ngβg)NTdV;]]>Fpg=∫ΩNpg(61ρgMg+1ρh∂sg∂shMh)RhdV+∫Ω▿Npg·KhkgρgμggdV+∫∂ΩNpgn·NqgqgdS;]]>CTT=∫ΩNT·(cTNT)dVdTdt;]]>KTT=∫Ω▿NT·(KT▿NT+(cwρwKhkwρwμw+cgρgKhkgρgμg)gNT)dV]]>KTTpw=∫Ω▿NT·(cwρwKhkwμw▿NpwpwNT)dVT;]]>KTTpg=∫Ω▿NT·(cgρgKhkgμg▿NpgpgNT)dVT;]]>FT=∫∂ΩNT·n·NqTqTdV+∫ΩNT·(-ΔHRh)dV;]]>式中：L为微分算子矩阵，D为弹性矩阵，ε为应变，εp为塑性应变，t为面力张量，f为体力张量，qT为热流，βh，βw和βg分别为水合物、水以及气体的热膨胀系数，Kh为水合物沉积物的固有渗透系数，kw和kg分别为水和气体的相对渗透系数，μw和μg分别是水和气体的粘滞性系数，Bw为水的体积膨胀系数，Mg、Mh以及Mw分别为气体、水合物和水的摩尔质量，Rh为水合物反应速率，δ为单位矩阵，n为法矢，g为重力加速度；联立上述方程，对指定的工况设置初值和边值条件，得到整个求解初值和边值问题的方程列式；要求解上述常微分方程组，首先进行时域离散，将常微分方程离散为代数方程系统来进行求解；常微分方程系统写成：0=L(U,U·,U··,t)-NF(U,t)Λ---(22)]]>0＝M(U,t)  (23)其中，L、NF、M分别代表刚度、约束力与约束，U为自由度；在求解该微分方程之前先化简拉格朗日乘子Λ；当0＝M(U,t)是线性的或时间相关时，且当约束力雅克比为常数时，消去约束M；其他情况下，则保留约束M进行求解；对于非线性问题，采用牛顿迭代法，其线性化公式表示为：EV··+DV·+KV=L-NFΛ---(24)]]>NV＝M  (25)其中，为刚度矩阵，称为阻尼矩阵，为质量矩阵，V为每一步计算得到的U的增量；假设自由度U通过如下关系化为一组无量纲化后的自由度表示：S是无量纲化系数矩阵；对线性方程组无量纲化后得到：其中，并且式中R是一个对角矩阵以使为1；假设U是某一时间步上的解，E是求解器在该时间步上对U的绝对误差做的估计，则表示为：(1MΣj1NjΣi(|Ei|Ai+R|Ui|)2)1/2<1---(28)]]>其中，Ai是自由度i的绝对误差，R为相对误差，M是参与耦合问题的物理场的个数，Nj是第j个物理场的自由度个数，该收敛准则等价于Ei<Ai+R|Ui|；(三)基于热力学方法的水合物沉积物临界状态本构模型的返回映射算法构造从弹性力学中应力σ与应变ε关系出发，计算试应力σtri和修正应力σcorr：σn+1=En+1:(ϵn+1-ϵn+1p)+Cn+1T=σtri+Δσcorr---(29)]]>其中，E为弹性模量；试应力为上一步传递过来的应力σn加上通过线弹性计算得到的应力增量：σtri=σn+En:Δϵn+1+CnΔTn+1+dEn+1:(ϵn-ϵnp)---(30)]]>修正应力是总的应力减去试应力得到的部分：Δσcorr=-En:Δϵn+1p+∂dE∂hΔh:(ϵn-ϵnp)+ΔCn+1Tn+1---(31)]]>其中，h是硬化参数矢量，下角标n+1表示第n+1步，C为热膨胀系数矢量，T为温度，E为弹性矩阵；弹性矩阵的变化由水合物对沉积物模量贡献给出：dEn+1=EnshχnΔsh+EnshΔsh∂χ∂hΔh---(32)]]>其中，χ是结构因子，sh为水合物饱和度；建立误差函数对误差函数进行离散：a=-ϵp+ϵnp+Δλr(σ,h,shn+1)=0---(33)]]>b＝‑h+hn+Δλv(σ，h，shn+1)＝0  (34)f(σ，h，shn+1)＝0  (35)其中，a,b分别为塑性应变矢量误差，f为屈服面函数，shn+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)f(σ(k),h(k),shn+1)+∂f∂σ(k):Δσ(k)+∂f∂h(k)·Δh(k)=0---(38)]]>公式(31)代入到(36)‑(38)中，得到：a(k)+En-1:Δσcorr(k)-En-1:((∂dE∂h·Δh(k)):(ϵn-ϵnp)+ΔCn+1Tn+1)+δλ(k)r(k)+Δλ(k)Δr(k)=0---(39)]]>b(k)‑Δh(k)+δλ(k)v(k)+Δλ(k)Δv(k)＝0  (40)f(k)+∂f∂σ(k):Δσ(k)+∂f∂h(k)·Δh(k)=0---(41)]]>Δr(k)=∂r∂σ(k):Δσ(k)+∂r∂h(k):Δh(k)---(42)]]>Δv(k)=∂v∂σ(k):Δσ(k)+∂v∂h(k):Δh(k)---(43)]]>对上述方程进行重新整理，得到：Δσ(k)Δh(k)=-A(k):a‾-δλ(k)A(k):r‾---(44)]]>其中：A-1(k)=En-1+Δλ(k)∂r∂σ(k)Δλ(k)∂r∂h(k)-En-1:∂dE∂h:(ϵn-ϵnp)Δλ(k)∂v∂σ(k)Δλ(k)∂v∂h(k)-I---(45)]]>a‾(k)=a(k)-En-1:ΔC(k)Tn+1b(k)---(46)]]>r‾(k)=r(k)v(k)---(47)]]>将应力和硬化参数增量代入f的表达式中解出δλ：-f(k)=∂f(k):Δσ(k)Δh(k)=-∂f(k):A(k):a‾(k)-δλ(k)∂f(k):A(k):r‾(k)---(48)]]>∂f(k)=[∂f∂σ(k),∂f∂h(k)]---(49)]]>δλ(k)=f(k)-∂f(k):A(k):a‾∂f(k):A(k):r‾(k)---(50)]]>模型中，屈服面f的表达式为：(p′-12γ(pcc′+pcd′+pcs′)R+R(1-γ2)pcc′)2((1-γ)p′+12γR(pcd′+pcs′)+pcc′)2+q2((1-α)Mp′+αMRρ′+RMpcc′)2=1---(51)]]>该模型引入应力比参数γ来控制干面和湿面的比例，由于受应力比的影响，屈服面将不再保持为椭圆，通过插值的方法定义转移应力ρ'；其中，M为临界状态线斜率，α是另一个与应力间距比相关的参数；p’为有效围压，q为广义剪应力，硬化参数h包括R，p′cd，p′cs，p′cc分别描述屈服前的非线性弹性，水合物胶接结构变化引起的剪胀行为，沉积物的固结过程以及胶接结构破坏引起的黏聚力降低；剪胀规律满足：D=dϵvp/dϵsp=(p′-12γ(pcc′+pcd′+pcs′)R+R(1-γ2)pcc′)·((1-α)Mp′+αMRρ′+RMpcc′)2q((1-γ)p′+12γR(pcd′+pcs′)+pcc′)2---(52).]]>