本發(fā)明屬于巖土工程領(lǐng)域，具體涉及一種考慮水合物分解的沉積物多場(chǎng)耦合模型的建模方法。

技術(shù)背景

天然氣水合物全球儲(chǔ)量巨大，分布廣泛，作為潛在的可替代能源，受到了世界各國(guó)政府的高度關(guān)注。然而，天然氣水合物大多賦存在永久凍土區(qū)和大陸邊緣，開(kāi)采難度大，其開(kāi)采過(guò)程會(huì)引起沉積層變形，不合理的開(kāi)采將導(dǎo)致地層沉降和海底滑坡，并可能誘發(fā)海嘯、溫室效應(yīng)等次生災(zāi)害。因此，天然氣水合物沉積層的物理力學(xué)特性研究顯得尤為重要。然而，水合物分解過(guò)程復(fù)雜，涉及多相多孔介質(zhì)內(nèi)含相變多場(chǎng)耦合問(wèn)題，目前國(guó)內(nèi)外缺少能夠精確描述水合物分解過(guò)程沉積層力學(xué)特性的多場(chǎng)耦合模型。

kimoto等首次嘗試將巖土力學(xué)與滲流、傳熱耦合，考慮了應(yīng)變率對(duì)變形計(jì)算的貢獻(xiàn)，分析了水合物分解過(guò)程對(duì)土體蠕變和強(qiáng)度衰減的影響，但是該模型極其復(fù)雜，參數(shù)眾多，且不能分析水合物賦存結(jié)構(gòu)變化引起的應(yīng)變軟化。rutqvist等基于tough+hydrate程序，利用flac商業(yè)軟件，進(jìn)行了滲流-力學(xué)半耦合計(jì)算，實(shí)現(xiàn)了巖土力學(xué)與水合物分解模擬的結(jié)合。該模型通過(guò)建立水合物飽和度與模量、強(qiáng)度的關(guān)系式，能夠模擬水合物分解對(duì)力學(xué)行為的影響，但不能模擬水合物分解引起的應(yīng)力松弛以及沉積層變形對(duì)傳熱和滲流的影響。klar等基于flac開(kāi)發(fā)了一套耦合傳熱、滲流和力學(xué)的程序，可以分析應(yīng)力松弛現(xiàn)象，但沒(méi)有考慮應(yīng)變軟化。國(guó)內(nèi)吳二林等、程遠(yuǎn)方等針對(duì)水合物降壓分解引起的工程問(wèn)題提出了相應(yīng)的多場(chǎng)耦合模型，但模型中都采用了相對(duì)簡(jiǎn)單的本構(gòu)關(guān)系，不能描述復(fù)雜力學(xué)行為。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的是為了克服上述現(xiàn)有模型的不足之處，提高對(duì)水合物分解過(guò)程沉積物力學(xué)特性的描述精度，提供一種考慮水合物分解的沉積物多場(chǎng)耦合模型及建模方法。本發(fā)明基于孔隙介質(zhì)力學(xué)基本理論，考慮了不同相物質(zhì)間的相互影響。從質(zhì)量守恒、動(dòng)量守恒和能量守恒出發(fā)，結(jié)合相關(guān)輔助方程，推導(dǎo)描述水合物分解過(guò)程沉積物物理力學(xué)行為變化的主要控制方程。利用comsolmultiphysics商用軟件的pde模塊以及結(jié)構(gòu)力學(xué)模塊中的外部材料接口，實(shí)現(xiàn)數(shù)值建模。其中，利用pde模塊分別建立液體連續(xù)性方程、氣體連續(xù)性方程、水合物分解控制方程以及能量守恒方程，利用結(jié)構(gòu)力學(xué)模塊中的外部材料接口，嵌入基于熱力學(xué)方法的水合物沉積物臨界狀態(tài)本構(gòu)模型。該模型參數(shù)皆可通過(guò)相應(yīng)實(shí)驗(yàn)獲得，便于實(shí)際應(yīng)用。

本發(fā)明的技術(shù)方案，通過(guò)以下步驟實(shí)現(xiàn)：

(1)基于孔隙介質(zhì)力學(xué)推導(dǎo)多場(chǎng)耦合控制方程；

(2)利用comsolmultiphysicspde模塊建模；

(3)采用基于熱力學(xué)方法的水合物沉積物臨界狀態(tài)本構(gòu)模型與返回映射算法；

(4)將本構(gòu)模型嵌入到comsolmultiphysics中實(shí)現(xiàn)水合物分解過(guò)程沉積物變形預(yù)測(cè)。

本發(fā)明的有益效果：模型基于孔隙介質(zhì)力學(xué)的基本概念，建立了溫度、孔隙壓力、應(yīng)力之間的耦合關(guān)系，理論基礎(chǔ)嚴(yán)謹(jǐn)，通用性強(qiáng)。利用comsolmultiphysics軟件中的pde模塊分別建立液體連續(xù)性方程、氣體連續(xù)性方程、水合物分解方程以及能量守恒方程，建模效率高，修改方便，利于擴(kuò)展至三維進(jìn)行實(shí)際工況模擬。利用結(jié)構(gòu)力學(xué)模塊中的外部材料接口，嵌入了基于熱力學(xué)方法的水合物沉積物臨界狀態(tài)本構(gòu)模型，使模型能夠描述水合物沉積物特有的力學(xué)性質(zhì)，如隨膠結(jié)結(jié)構(gòu)破壞而引起的應(yīng)力軟化、剪脹、屈服面上剪脹部分與簡(jiǎn)縮部分的不等性以及屈服面形狀等。同時(shí)，采用向后歐拉積分形式的返回映射算法，利用c++語(yǔ)言編寫模型程序，算法通用性強(qiáng)，程序可擴(kuò)展性好。

附圖說(shuō)明

圖1是本發(fā)明的考慮水合物分解的沉積物多場(chǎng)耦合模型原理圖。

圖2是本發(fā)明的返回映射算法流程圖

圖3(a)是本構(gòu)模型模擬應(yīng)力應(yīng)變曲線圖。

圖3(b)是本構(gòu)模型模擬剪脹曲線圖。

圖4是模型驗(yàn)證采用的網(wǎng)格和尺寸圖。

圖5是孔隙壓力模擬結(jié)果圖。

圖6是水合物飽和度模擬結(jié)果圖。

圖7是溫度模擬結(jié)果圖。

圖8是有效應(yīng)力路徑模擬結(jié)果圖。

具體實(shí)施方式

以下結(jié)合附圖和技術(shù)方案，進(jìn)一步說(shuō)明本發(fā)明的具體實(shí)施方式。

實(shí)施例

一種考慮水合物分解的沉積物多場(chǎng)耦合模型的建模方法，步驟如下：

(1)基于孔隙介質(zhì)力學(xué)推導(dǎo)多場(chǎng)耦合控制方程

根據(jù)孔隙介質(zhì)力學(xué)，α相物質(zhì)體積分?jǐn)?shù)為

其中，α為s、w、g或h，分別代表沉積物顆粒、水、氣體和水合物，v表示體積，vα表示α相物質(zhì)的體積；水合物、孔隙水和孔隙氣飽和度分別定義為sh＝nh/n、sw＝nw/n和sg＝ng/n，其中，n是孔隙率。三者之間關(guān)系如下所示：

sw+sh+sg＝1(2)

根據(jù)孔隙介質(zhì)力學(xué)，α相物質(zhì)的質(zhì)量守恒方程表示為：

其中，為速度散度，ρα表示α相物質(zhì)的密度；dmα/dt表示α相物質(zhì)的質(zhì)量累積，即源或匯。

對(duì)于沉積物顆粒，開(kāi)采過(guò)程中的溫度和孔壓變化不足以產(chǎn)生較大的沉積物顆粒膨脹或收縮，因此這里忽略溫度和孔壓對(duì)沉積物顆粒密度的影響。沉積物顆粒的連續(xù)性方程為：

以沉積物顆粒作為骨架，則有εv為體積應(yīng)變；

對(duì)于水和氣體，其連續(xù)性方程為：

水合物連續(xù)性方程為：

沉積物骨架由沉積物顆粒與水合物共同構(gòu)成，且不考慮水合物從沉積物骨架中擠出，故vs＝vh；

方程中的質(zhì)量累積項(xiàng)滿足化學(xué)反應(yīng)過(guò)程中的質(zhì)量守恒，因此

將流速使用達(dá)西流速表示vw＝qw/nw+vs和vg＝qg/ng+vs分別代入到連續(xù)性方程(6)和(7)中，其中，vw、vg和vs分別是孔隙水、孔隙氣以及土骨架的速度，并且忽略nw、ng和密度的空間分布對(duì)質(zhì)量守恒的影響，得：

通過(guò)孔隙介質(zhì)力學(xué)理論描述水合物沉積物的能量守恒；為了簡(jiǎn)化方程，認(rèn)為外界對(duì)物體做功主要轉(zhuǎn)化為動(dòng)能，滿足機(jī)械能守恒定律；同時(shí)，為了提高計(jì)算速度，不考慮內(nèi)力做功對(duì)溫度變化的微弱影響。方程如下所示：

其中，導(dǎo)熱系數(shù)kt利用疊加原理獲得：kt＝nskts+ngktg+nhkth+nwktw；kts、ktg、kth和ktw分別為固體、氣體、水合物和水的導(dǎo)熱系數(shù)；ct是水合物沉積物的總熱容，滿足ct＝ρsnscs+ρgngcg+ρhnhch+ρwnwcw，其中，cs、cg、ch、cw分別為沉積物顆粒、氣體(需壓力修正)、水合物和水的比熱；δh為每mol水合物相變引起的熱量變化，通過(guò)kamath回歸方程計(jì)算，滿足δh＝c1+c2t，c1和c2為兩個(gè)回歸參數(shù)。

力平衡方程的具體描述為：

其中，σ為總應(yīng)力，fα為α相物質(zhì)流體的滲透阻力；水合物沉積物為多孔介質(zhì)，有效應(yīng)力是產(chǎn)生變形和影響強(qiáng)度的主要因素；根據(jù)有效應(yīng)力原理，截面上應(yīng)力等于有效應(yīng)力與孔隙壓力之和，表示為：

σ＝σ'-ppδ(16)

其中，σ'為有效應(yīng)力，pp為孔隙壓力；δ為單位矩陣；對(duì)于水合物沉積物，其孔隙壓力表達(dá)式如下所示：

其中，pw和pg分別為孔隙水壓力和孔隙氣壓力。

(2)comsolmultiphysicspde建模

控制方程包括方程(8)、(10)-(13)，其中偏微分方程通過(guò)有限元離散化為常微分方程進(jìn)行求解，采用的單元類型為混階lagrange單元；選取lagrange插值形函數(shù)nu，npw，npg，nt，nqw，nqg，nqt代入到控制方程中，并通過(guò)變分原理得到常微分方程：

kuu:u+kupgpg+kupwpw＝fu(18)

其中，

式中：l為微分算子矩陣，d為彈性矩陣，ε為應(yīng)變，ε^p為塑性應(yīng)變，t為面力張量，f為體力張量，qt為熱流，βh，βw和βg分別為水合物、水以及氣體的熱膨脹系數(shù)，kh為水合物沉積物的固有滲透系數(shù)，kw和kg分別為水和氣體的相對(duì)滲透系數(shù)，μw和μg分別是水和氣體的粘滯性系數(shù)，bw為水的體積膨脹系數(shù)，mg，mh以及mw分別為氣體、水合物和水的摩爾質(zhì)量，rh為水合物反應(yīng)速率，δ為單位矩陣，n為法矢，g為重力加速度。

聯(lián)立上述方程，對(duì)指定的工況設(shè)置初值和邊值條件，得到整個(gè)求解初值和邊值問(wèn)題的方程列式。

要求解上述常微分方程組，首先進(jìn)行時(shí)域離散，將常微分方程離散為代數(shù)方程系統(tǒng)來(lái)進(jìn)行求解。通常常微分方程系統(tǒng)寫成：

0＝m(u,t)(23)

其中，l、nf以及m分別代表剛度、約束力與約束，u為自由度；在求解該微分方程之前先化簡(jiǎn)拉格朗日乘子λ；當(dāng)0＝m(u,t)是線性的或時(shí)間相關(guān)時(shí)，且當(dāng)約束力雅克比為常數(shù)時(shí)，算法消去約束m，其他情況下，則保留約束m進(jìn)行求解。

對(duì)于非線性問(wèn)題，采用牛頓迭代法，其線性化公式表示為：

nv＝m(25)

其中，為剛度矩陣，稱為阻尼矩陣，為質(zhì)量矩陣，v為每一步計(jì)算得到的u的增量。

由于不同物理場(chǎng)控制方程之間剛度差別較大，故為了降低條件數(shù)，減小剛度的奇異性，需要進(jìn)行無(wú)量綱化處理。假設(shè)自由度u可以通過(guò)如下關(guān)系化為一組無(wú)量綱化后的自由度表示：

s是無(wú)量綱化系數(shù)矩陣；對(duì)線性方程組無(wú)量綱化后得到：

其中，并且式中r是一個(gè)對(duì)角矩陣以使的大小為1。

對(duì)于耦合氣-水兩相流，由于氣體滲透性和壓縮性與水的滲透性和壓縮性相差甚大，如果以氣壓和水壓作為獨(dú)立未知量會(huì)使得剛度陣奇異，求解困難增加。故這里采用水壓和有效水飽和度作為未知量，并在剛度矩陣中引入松弛因子以避免剛度矩陣的奇異，有效且穩(wěn)定地求解壓力。

comsol使用相對(duì)誤差控制積分精度，利用絕對(duì)誤差控制誤差累積；假設(shè)u是某一時(shí)間步上的解，e是求解器在該時(shí)間步上對(duì)u的絕對(duì)誤差做的估計(jì)，那么當(dāng)誤差滿足下面條件就認(rèn)為計(jì)算在該時(shí)間步上收斂了：

其中，ai是自由度i的絕對(duì)誤差，r為相對(duì)誤差，m是參與耦合問(wèn)題的物理場(chǎng)的個(gè)數(shù)，nj是第j個(gè)物理場(chǎng)的自由度個(gè)數(shù)，該收斂準(zhǔn)則等價(jià)于ei<ai+r|ui|。

(3)基于熱力學(xué)方法的水合物沉積物臨界狀態(tài)本構(gòu)模型的返回映射算法構(gòu)造

采用向后歐拉法建立應(yīng)力計(jì)算的返回映射算法，包括兩部分：第一部分是將本構(gòu)方程離散成代數(shù)方程的積分算法，第二部分是求解非線性代數(shù)方程的迭代算法；具體實(shí)現(xiàn)步驟包括兩步：首先通過(guò)彈性預(yù)測(cè)計(jì)算一個(gè)試應(yīng)力，理論上該試應(yīng)力偏離屈服面，然后通過(guò)塑性調(diào)整將應(yīng)力拉回到屈服面上。

從彈性力學(xué)中應(yīng)力的定義出發(fā)，計(jì)算出試應(yīng)力和修正應(yīng)力：

其中，試應(yīng)力為上一步傳遞過(guò)來(lái)的應(yīng)力加上通過(guò)線彈性計(jì)算得到的應(yīng)力增量：

修正應(yīng)力是總的應(yīng)力減去試應(yīng)力得到的部分：

其中，h是硬化參數(shù)矢量，下角標(biāo)n+1表示第n+1步，c為熱膨脹系數(shù)矢量，t為溫度，e為彈性矩陣。彈性矩陣的變化由水合物對(duì)沉積物模量貢獻(xiàn)給出：

其中，χ是結(jié)構(gòu)因子。

建立誤差函數(shù)對(duì)誤差函數(shù)進(jìn)行離散：

b＝-h+hn+δλv(σ，h，shn+1)＝0(34)

f(σ，h，shn+1)＝0(35)

其中，a,b分別為塑性應(yīng)變矢量誤差，f為屈服面函數(shù)，shn+1為第n+1步水合物飽和度，因?yàn)榉e分格式為隱式，所以第n+1步的水合物飽和度直接由質(zhì)量守恒方程獲得；r為塑性應(yīng)變?cè)隽康姆较颍瑅為硬化參數(shù)矢量的方向。

通過(guò)離散得到線性代數(shù)方程，獲得第k步的線性代數(shù)方程：

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)

對(duì)上述方程進(jìn)行重新整理，得到：

其中：

將應(yīng)力和硬化參數(shù)增量代入f的表達(dá)式中解出δλ：

模型中，屈服面f的表達(dá)式為：

該模型引入應(yīng)力比參數(shù)γ來(lái)控制干面和濕面的比例，由于受應(yīng)力比的影響，屈服面將不再保持為橢圓，通過(guò)插值的方法定義轉(zhuǎn)移應(yīng)力ρ'。其中，m為臨界狀態(tài)線斜率，α是另一個(gè)與應(yīng)力間距比相關(guān)的參數(shù)。

p’為有效圍壓，q為廣義剪應(yīng)力，硬化參數(shù)h包括r，p′cd，p′cs，p′cc分別描述屈服前的非線性彈性，水合物膠接結(jié)構(gòu)變化引起的剪脹行為，沉積物的固結(jié)過(guò)程以及膠接結(jié)構(gòu)破壞引起的黏聚力降低。

剪脹規(guī)律滿足：

圖2為具體的算法流程圖。圖3為模型預(yù)測(cè)結(jié)果與三軸試驗(yàn)結(jié)果的對(duì)比，結(jié)果表明本模型對(duì)試驗(yàn)的應(yīng)力應(yīng)變關(guān)系和剪脹關(guān)系有較好的擬合效果。

(4)將本構(gòu)模型嵌入到comsolmultiphysics中實(shí)現(xiàn)水合物分解過(guò)程沉積物變形預(yù)測(cè)。

comsolmutiphysics結(jié)構(gòu)力學(xué)模塊提供了可以進(jìn)行本構(gòu)二次開(kāi)發(fā)的接口。該接口共設(shè)置了三類變量進(jìn)行數(shù)據(jù)傳遞，分別為輸入變量、輸出變量與狀態(tài)變量。本模型中所涉及的輸入變量為應(yīng)變和溫度，輸出變量為應(yīng)力，狀態(tài)變量包括r，p′cd，p′cs，p′cc以及塑性應(yīng)變，材料參數(shù)包括基于熱力學(xué)方法的水合物沉積物臨界狀態(tài)本構(gòu)模型參數(shù)m，α，γ以及水合物飽和度sh。將接口對(duì)應(yīng)的變量和參數(shù)傳遞給本構(gòu)模型c++程序中，實(shí)現(xiàn)變形計(jì)算。

通過(guò)驗(yàn)算uchida等人論文中的標(biāo)準(zhǔn)考題對(duì)模型進(jìn)行驗(yàn)證。算例中初始溫度設(shè)置為284k，對(duì)應(yīng)相平衡曲線上的壓力值為8.3mpa。幾何尺寸、初始條件、邊界條件如圖4所示，考察兩個(gè)局部位置r/rwell＝1.5和5.0處的孔壓pp、水合物飽和度sh、溫度t以及相關(guān)力學(xué)參數(shù)的變化。初始孔隙壓力、豎向有效應(yīng)力和水平有效應(yīng)力分別設(shè)置為13mpa、3mpa和1.5mpa，對(duì)應(yīng)的初始應(yīng)力比為k0＝0.5。水合物初始飽和度為0.5，氣體飽和度為殘余氣體飽和度，壓力從13mpa經(jīng)過(guò)兩天降低到7mpa，并在之后的18天中維持7mpa不變。圖5-圖8為本模型與uchida等人使用flac計(jì)算得到的標(biāo)準(zhǔn)考題解的比較，可以發(fā)現(xiàn)二者基本吻合。圖5-圖7表明：距離開(kāi)采井較近的區(qū)域孔壓和溫度降低速度較遠(yuǎn)處更快，水合物分解量更大。圖8為有效應(yīng)力路徑的計(jì)算結(jié)果。由于初始條件為k0固結(jié)，因此在降壓初始階段，沉積層從k0固結(jié)向各向同性過(guò)渡，導(dǎo)致偏應(yīng)力降低；同時(shí)，由于上下約束，體積收縮引起有效小主應(yīng)力降低，而有效大主應(yīng)力隨降壓而升高，導(dǎo)致偏應(yīng)力再次升高。水合物分解會(huì)引起應(yīng)力松弛，不論是平均有效應(yīng)力還是偏應(yīng)力都會(huì)降低。