天然氣水合物試采中節流螺旋管段微米級砂粒運移沉積規律數值模擬

蘇義腦 陳 燁 孫曉峰 閆 鐵 曲晶瑀 段瑞溪

1.東北石油大學石油工程學院 2.中國石油海洋工程有限公司

0 引言

我國海域天然氣水合物（以下簡稱水合物）資源豐富，其儲層特點與日本類似，儲層埋深淺，膠結較弱，未固結成巖，水合物在砂土中以填充、骨架和膠結幾種土體賦存方式存在[1-5]，一旦試采破壞了水合物的賦存條件[6-7]，起膠結和骨架支撐作用的水合物土體就會失去支撐，儲層砂粒在地層流體的攜帶下流向井筒。水合物儲層巖性以微米級黏土、粉砂、細砂為主，南海多個站位巖心樣品粒度分布測試結果顯示粒徑小于63 μm的粉砂超過65%[8-10]，有的站位超過70%[11]，甚至有測試結果顯示小于63 μm的黏土和粉砂占比接近94%[12]。對于微米級的泥質粉砂，防砂技術難度較大，目前水合物試采所用的防砂技術主要包括機械防砂、加裝防砂網、射孔砂篩防砂、Geoform防砂系統等，這些防砂技術效果并不理想，實踐和研究表明，防砂系統很難對粒徑44 μm以下的砂粒發揮作用[13]。近年來，陸地（加拿大，2007年）和海域（日本，2013年、2017年）水合物開采，都發生了因微米級砂粒突破防砂設施進入井筒，引起井筒砂埋、堵塞，迫使水合物試采作業提前終止的案例[14-16]。日本共進行了兩次大規模海域水合物開采試驗，其中2013年為世界首次海域水合物開采試驗，共產氣6天，日產氣2×104m3，后由于大量儲層砂粒突破防砂設施進入井筒，導致試采提前結束[15,17-19]。2017年，日本總結第一次水合物試采經驗，進行第二次海域水合物試采作業，同樣由于井筒沉砂堵塞而終止產氣[20]。因而，水合物開采過程中不僅需要考慮儲層微米級砂粒突破井筒防砂設施后井筒設備的磨損問題，更要考慮微米級砂粒在井筒內的運移、沉積、堵塞問題，需要對微米級砂粒進入井筒后的運移沉積規律特征進行研究。

目前，水合物出砂研究主要集中在地質儲層中，包括儲層巖石力學特性（產砂機理）[21-27]、驅動力[28]、出砂風險[29-33]、評價防砂方法有效性[34-37]幾個方面，而微米級砂粒隨地層流體進入井筒后的運移、沉積、堵塞研究尚未見文獻報道。本文研究的目的是：①通過數值模擬，獲得微米級砂粒在節流管線螺旋段的運移沉積特征；②獲得粒徑、水速等參數對微米級砂粒運移沉降的影響規律；③獲得不同情況下的砂粒臨界不沉積水速，建立局部復雜井段微米級砂粒運移沉積預測模型，為水合物開采井筒砂埋提供預測方法；④為減小砂堵風險、確定合理開采工作制度、保證井筒安全提供依據，為后續水合物儲層適度防砂工藝的發展提供技術支撐。

1 計算方法

1.1 計算模型

降壓開采法是目前水合物開采的主流方法之一。前蘇聯麥索亞哈、加拿大麥肯齊、美國阿拉斯加北坡、日本南海海槽東部、日本南海海槽、中國南海神狐都采用過降壓法開采水合物[20]，日本2013年在南海海槽東部采用降壓開采法進行了第一次海域水合物試采作業，其開采系統如圖1-a所示。

降壓破壞了水合物穩定賦存的壓力條件，使水合物發生分解，產生甲烷氣和水，分解后產生的氣液兩相在井下分離并通過獨立的管路產出，其中的甲烷等烴類氣體通過油管通道舉升到地面，分離水通過節流管線單獨輸運到井口。分離水在井下分離后經過C型通道進入節流管線，由于節流管線下部為剛性管路，且上下接箍位置不在同一側，因而存在一段螺旋段，如圖1-b所示。可見，水流通道既存在C型段，也存在螺旋段，是流道復雜砂粒易沉降的位置。為此，對水合物開采局部復雜管路微米級砂粒運移規律進行數值模擬研究，建立計算的幾何模型如圖1-c所示。由于節流管線為水流通道，為此選定水為連續相，微米級砂粒為分散相。

1.2 控制方程

連續性方程為：

式中下標l和s分別表示水和微米級砂粒；α表示濃度；表示速度向量；ρ表示密度。

每個計算單元相互滲透，因此他們濃度的和是一致的。水和砂粒的動量方程可以表示為：

圖1 水合物開采系統、節流管線及計算模型示意圖

當αl≤0.8時，水—砂粒的交換系數可以表示為：

當αl＞0.8時，水—砂粒的交換系數可以表示為：

拖拽系數（CD）可表示為：

砂粒的雷諾數定義為：

式中ds表示砂粒粒徑。

各相的應力—應變張量可表示為：

Realizablek—ε模型中的湍動能（k），耗散率（ε）為：

式中σk和σ?分別表示湍動能方程k和耗散率方程ε的紊流普朗特數；C2表示常數。

Gk表示由于層流速度梯度而產生的湍流動能，可以表示如下：

1.3 邊界條件及網格劃分

水和均勻分散的微米級砂粒從節流管線底部進入，經C型段、螺旋段后進入垂直段而后到達地面，如圖2所示。在入口處選擇速度入口邊界條件，出口選擇壓力出口邊界條件，壁面采用無滑移壁面條件。

適當的網格劃分軟件和高精度網格質量是CFD計算的必要先決條件。為提高計算準確性，加快計算和收斂速度，使用ICEM軟件采用“自下而上”方式對計算流場劃分六面體結構網格，圖2為整體計算流場及局部位置三維網格劃分示意圖。設置邊界層網格，管路徑向第一層網格高度的選取符合y＋≈30，網格增長率取1.2，沿徑向第一層網格高度采用式（16）計算

式中Δy表示第一層網格高度，mm；Re表示雷諾數；D表示水力特征長度，mm。

為了檢測網格獨立性，劃分3種數量網格（158 347，200 816，411 075）用于流場（44 μm，1%濃度）模擬，模擬結果如表1所示。可以看出其相對誤差均小于4.5%。綜合考慮計算資源消耗、計算準確性，選擇網格數200816進行數值模擬計算。

通過有限體積法對控制方程進行離散化之后，將相位耦合的SIMPLE算法應用于壓力—速度耦合。動量方程使用二階迎風格式進行離散化以提高計算精度。對每個尺度殘差分量采用10－5的收斂準則來確定兩個連續迭代之間的相對誤差。時間步長取0.005 s，每個步長迭代20步。

選取3種粒徑砂粒（44 μm，10 μm，4 μm），3種出砂濃度（1%，5%，10%），多種流速進行數值模擬計算，獲得不同流速下微米級砂粒在復雜井段的運移沉積情況，以進一步確定不同粒徑砂粒在不同濃度下運移的臨界不沉積水速。

圖2 整體計算區域及局部位置3D網格劃分圖

表1 網格獨立性檢驗表

2 模擬結果及分析

2.1 微米級砂粒沉積分析

圖3顯示了44 μm粒徑、10%出砂濃度在不同水速下微米級砂粒濃度分布云圖。

由圖3可以看出水速較小時，微米級砂粒在管路中的沉積較為明顯，主要集中在C型段和螺旋段的下部（圖3-a）；局部砂粒濃度最高達0.55。隨著流速的增加，砂粒沉積情況改善較大，C型段、螺旋段下部沉積減少，局部砂粒濃度最大為0.24，集中在螺旋段中上部，管路其余部分多為0.12左右（圖3-b），繼續增大水速，微米級砂粒沉積越來越小，只在C型段轉角留存有部分沉積（圖3-c）。可見，微米級砂粒對水速較為敏感，較小的水速變化就能引起較大的砂粒濃度的變化。微米級砂粒較易沉積的部位主要在C型段的拐角處及螺旋段，這是由于流向變化和重力作用引起的。

圖3 44 μm粒徑、10%出砂濃度在不同水速下微米級砂粒濃度分布云圖

2.2 微米級砂粒臨界不沉積流速的確定

對3種粒徑、3種出砂濃度下微米級砂粒在復雜井段內的運移沉積情況進行數值模擬，共獲得100多組模擬數據，模擬結果如圖4所示。

圖4 3種粒徑、3種濃度下流速與復雜管段內砂粒沉積濃度關系圖

以粒徑44 μm顆粒1%濃度下，管路內微米級砂粒濃度隨水速變化規律為例，由圖4-a可以看出，當水速較大時，復雜管段內砂粒濃度無限趨近于地層出砂濃度1%，說明管路內幾乎沒有砂粒沉降，曲線變化平緩。隨著水速的逐漸減小，管路內微米級砂粒濃度逐漸增加，當水速減小到一定程度后，小幅度的速度減小將引起管路內砂粒濃度的急劇升高，濃度呈指數增長趨勢。因此可將曲線變化大體劃分為3個部分：平緩區、指數增長區以及兩者交匯的過渡區。臨界不沉積水速位于兩種變化趨勢交匯的過渡區內。如何確定臨界不沉積水速，筆者提出了一種方法：由于臨界不沉積水速位于過渡區內，可分別對平緩區和指數增長區進行直線擬合，擬合曲線分別表示兩種區域的變化趨勢，則這兩條直線的交點即為臨界點，其橫坐標即為臨界不沉積水速。水速大于這個值，管路內砂粒沉積量較小，水速小于這個值，管路內的砂粒沉積量將急劇增加，且水速越小，沉積量越大。由此可得到44 μm，1%濃度下微米級砂粒在水合物試采局部復雜井段的臨界不沉積水速為0.106 5 m/s。同理，可得到其他工況下的臨界不沉積水速，如圖4-b～i所示。

2.3 臨界不沉積流速結果分析

微米級砂粒臨界不沉積水速隨粒徑及地層出砂濃度的變化規律如圖5、6所示。

圖5 不同地層出砂濃度下臨界流速隨粒徑變化規律圖

由圖5可以看出，隨著粒徑的增大，微米級砂粒的臨界不沉積水速逐漸增大，粒徑較小時，出砂濃度的變化對臨界不沉積水速的影響較小，粒徑較大時，出砂濃度差別造成的臨界不沉積水速的變化變得更加明顯。

圖6 不同粒徑下臨界流速隨地層出砂濃度變化規律圖

由圖6可以看出，隨著地層出砂濃度的增加，微米級砂粒的臨界不沉積水速逐漸增加，且濃度較大時，粒徑變化造成的臨界不沉積水速的變化更大。造成這種現象的原因是因為隨著粒徑的增大，自重增加，砂粒懸浮需要更大的流速，另一方面，砂粒較小時，比表面積更大，砂粒更易懸浮在溶液中，不易沉降。而且地層出砂濃度越大，井筒內含砂量越多，沉積更容易發生。粒徑和地層出砂濃度的雙重作用導致粒徑較大時，濃度差別造成的臨界不沉積水速變化更加明顯，出砂濃度較大時，粒徑差別造成的臨界不沉積水速變化更大。

2.4 微米級砂粒沉積預測模型建立

由于試采局部井段管路復雜，室內實驗存在較大難度，數值模擬計算需要耗費較長的時間和計算資源，為獲得一種相對便捷且成本較低的計算方式，對上述模擬數據進行分析，水合物試采局部復雜管路微米級砂粒的濃度（C）是地層出砂濃度（C0）、液流速度（v）、液流密度（ρl）、液流黏度（μl）、砂粒粒徑（ds）及螺旋段井筒直徑（Dlx）的函數，對這些變量進行分析，如表2所示，這些變量之間的數學關系可表示為：

表2 獨立變量分析表

應用Buckingham-Π定理，對上述變量進行無量綱分析，組合成無量綱群，將復雜問題簡化為較小組合變量之間的關系。為此定義4個具有實際意義的無量綱量，對所研究變量進行無量綱化。第1個無量綱量為復雜管路井筒內微米級砂粒濃度（π1，C），是因變量，其他變量為自變量；第2個無量綱量為地層出砂濃度（π2，C0），表示地層出砂速率對砂粒沉積運移的影響；第3個無量綱量為（π3，ds/Dlx），代表了粒徑對微米級砂粒運移沉積的影響；第4個無量綱量為（π4，ρ1vDlx/μ1）代表了螺旋管雷諾數，式中的特征長度為螺旋管的水力直徑，具體如表3所示。

表3 獨立變量無量綱化表

4組具有實際物理意義的無量綱量之間的數學關系可以描述如下：

利用OriginPro 2019軟件非線性擬合工具對其進行分析，確定擬合方程如下：

式中a1、a2、a3、a4表示系數。對數據進行擬合，擬合結果如表4所示。

則最終的復雜管路微米級砂粒沉積預測方程為：

對數據點模型預測值進行統計誤差分析來檢驗模型的準確性。檢驗模型準確性的統計學參數包括：平均誤差（ME）、平均絕對誤差（MAE）、標準差（SD）、平均絕對百分比誤差（MAPE）、均方根誤差（RMSE），這些統計參數定義如式（21）～（25）所示[39]。統計參數分析結果如表5所示。可以看出ME、MAE、RMSE均在0.7以下，SD和MAPE也在9.1以內，結合表4中的相關系數、殘差平方和及擬合優度，可以看出模型誤差較小。

表5 統計參數分析結果表

式中n表示模擬的組數；ecal、esim分別表示計算、數值模型的井筒砂粒濃度。

x的表達式為：

為了更好地描述水合物試采局部復雜井段沉砂的嚴重程度，方便計算和預測不同條件下的砂粒臨界不沉積水速，定義砂沉積濃度比這一概念，即：試采后井筒內的含砂濃度C與地層出砂濃度C0（與降壓方式、地層出砂速率等有關）之比。即

當e＝1時，說明進入井筒內的砂量等于運移出井筒的砂量，此時的水速即對應為臨界不沉積水速；當e＞1時，說明井筒內有沉積，且e值越大，沉積越嚴重。

定義砂沉積濃度比后，利用式（20），可計算得到井筒內的沉砂濃度，進一步利用式（27）計算砂沉積濃度比，直觀顯示井筒內的沉砂情況。另一方面，令e＝1，可以反推計算得到不同情況下復雜管路中的砂粒臨界不沉積水速，從而為安全合理安排水合物生產制度、降壓幅度提供參考。

3 結論

1）螺旋段內微米級砂粒沉積情況隨著水速的增加而逐漸改善，其中螺旋段上部的砂粒清潔難度要大于螺旋段下部。

2）提出了一種微米級砂粒臨界不沉積水速的確定方法，得到了3種粒徑3種濃度下水合物試采局部復雜管路微米級砂粒的臨界不沉積水速。

3）微米級砂粒臨界不沉積水速隨著粒徑和地層出砂濃度的增大而增大，且粒徑較大時濃度差別造成的臨界不沉積水速變化更加明顯，濃度較大時粒徑差別造成的臨界不沉積水速變化更大。

4） 提出了砂沉積濃度比的概念，當濃度比為1時，井筒內沒有砂粒沉積，當濃度比大于1時，井筒內有沉積，且濃度比越大沉積越嚴重。

5）得到了水合物試采局部復雜井段井筒沉砂濃度預測模型，計算不同條件下的砂粒臨界不沉積水速，為安全合理安排水合物生產制度、降壓幅度提供有益參考。