混合無機納米孔內流體運移規律研究*

高明玉，劉丹，董志立，梁爽**，王京坤，羅鑫

(1. 東北石油大學石油工程學院，黑龍江大慶 163318；2. 東北石油大學提高油氣采收率教育部重點實驗室，黑龍江大慶 163318；3. 國家管網集團油氣調控中心，北京 100013；4. 大慶油田第八采油廠，黑龍江大慶 163514；5. 大慶油田第四采油廠，黑龍江大慶 163511)

隨著社會經濟發展，對油氣資源的需求日益增大，非常規油氣資源越來越受到重視，頁巖油也隨之成為世界各國油氣勘探的焦點。早在2012 年，美國便依靠大力開發頁巖油實現了原油產量劇增，2019 年美國超過俄羅斯和沙特阿拉伯成為全球最大的原油生產國。由此可見，頁巖油的勘探開發對世界能源格局有著不可忽視的影響[1]。

頁巖油是以頁巖為主的頁巖層系中所含的原地滯留油氣資源，圈閉界限不明顯，無法形成自然工業產能[2]。與常規儲層不同，頁巖儲層具有極低的滲透率和豐富的納米孔隙，對于幾納米到幾十納米的微小孔隙，傳統流體力學理論和試驗方法難以描述流體的吸附、運移等細節，基于牛頓力學的分子動力學模擬則能很好地解決這些問題[3]。分子動力學模擬以分子為研究對象，將整個系統作為有一定特征的分子集合，利用牛頓運動定律對系統內分子的運動軌跡進行模擬，并通過系綜計算體系的結構和性質。盡管頁巖油氣賦存及運移規律的研究已經取得了諸多成果，但頁巖儲層礦物成分復雜，存在大量多種成分混合的納米孔，對于混合孔隙內流體運移規律的相關研究較少，仍需進一步分析。

筆者以辛烷和水為研究對象，采用分子動力學模擬方法開展了石英和伊利石混合的無機質孔隙內流體賦存狀態和運移規律的研究，通過分析孔隙內流體的密度分布、速度分布、滑移長度、黏度等參數，分別討論了孔徑大小、溫度、壓力梯度及含水率對納米孔隙內流體賦存和運移的影響。

1 模型構建及模擬方案

1.1 模型構建

頁巖基質骨架由有機質和無機質組成，有機質主要成分為干酪根，無機質主要成分為石英和黏土礦物，黏土礦物中以伊利石、綠泥石及伊蒙混層為主要成分，部分含有極少量綠蒙混層，其中伊利石占比最大，平均含量達到74.1%[4]，因此，以石英與伊利石組合模型構建頁巖孔隙。伊利石是以Al原子為中心的八面體層，夾在2 個以Si 原子為中心的四面體之間作為2∶1 型黏土礦物。在每個復制單元晶胞中，一個Si4+被一個Al3+取代，形成同構取代。層間反陽離子(鉀離子，K+)在伊利石層間空間中隨機分布，以抵消伊利石層內同構取代引起的靜電電荷，伊利石結構原子坐標及晶胞參數分別見表1、表2。石英是一種由SiO2四面體通過共同頂點連接成的三維骨架結構。在石英晶體中，硅原子位于四面體中心，1 個硅原子與4 個氧原子成鍵，同時1 個氧原子與2 個硅原子成鍵，氧原子處于硅氧四面體的公用頂點上。

表1 伊利石模型原子坐標

表2 伊利石模型晶胞參數

晶體的形貌是晶體晶面生長發育的結果，在晶體生長發育的過程中，生長速度快的晶面往往消失，生長速度慢的晶面容易保留下來，成為晶體形貌最主要的晶面。在不破壞Si—O 鍵的情況下，伊利石的(001)面是伊利石形貌最重要的晶面，該晶面用于考察外界與晶面的相互作用，有足夠的代表性。Chilukoti 等[5]通過試驗發現液態烷烴更傾向于吸附在α-SiO2的(100)晶面。因此，模型構建將α-SiO2的(100)晶面作為SiO2的代表面，與伊利石的(001)面組合成伊利石-二氧化硅模型，模型結構見圖1；狹縫孔隙由兩個平行的伊利石-二氧化硅組合壁面構成，向兩個壁面中添加一定數量的n-C8H18[6-7]，構成孔徑分別為2，5，8 nm 寬的孔隙。模型在3個方向均采用周期性邊界條件，為了消除邊界效應，在3 個方向上分別添加了0.1 nm 的真空區域，伊利石-二氧化硅混合頁巖孔隙內流體分布初始結構見圖2。

圖1 伊利石-二氧化硅組合模型結構

圖2 伊利石-二氧化硅混合頁巖孔隙內流體分布初始結構

1.2 模擬方案

為了保證模擬結果的準確性，分子動力學模擬過程中選取的勢能參數至關重要。采用SPC/E 模擬水分子，ClayFF 力場模擬石英和伊利石基質[7-8]，此力場廣泛應用于黏土界面模擬。ClayFF力場僅考慮水分子、羥基、可溶性多原子分子和離子間的鍵伸縮和鍵角彎曲項，其他所有的相互作用力均采用非鍵相互作用勢即Lennard-Jones勢與庫侖力之和進行模擬?？紫秲鹊膎-C8H18分子采用OPLS-AA(Optimized Potentials for Liquid Simulation)全原子力場模擬[7]，其特點是對于不同原子間的勢能參數通過試驗結果擬合得到，此力場被廣泛應用于烷烴、聚合物、生物大分子等，具有很高的可靠性。短程非鍵范德華相互作用使用Lennard-Jones(12-6)表示，其截斷半徑設置為1.2 nm，符合截斷半徑不能大于最小模擬盒子尺寸一半的規則。采用基于傅里葉的Ewald求和方法計算了長程靜電力相互作用，粒子與粒子間的粒子網格采用PPPM方法，精度值設置為10-6。通過Parrinello-Rahman穩壓器控制體系的壓力，通過Nose Hoover恒溫器控制溫度，模擬初始在等壓等溫NPT系綜下進行，溫度360 K，壓力30 MPa，體積波動僅設置在z方向。設置模擬的時間步長為1 fs，在NPT系綜下弛豫1000 ps，取最后500 ps數據用來分析。

2 結果分析

2.1 孔喉尺寸對流體運移影響

溫度為360 K，孔徑分別為2，5，8 nm 狹縫內流體的連續密度分布曲線及壓力梯度為0.002 kcal/(mol·?)(換算單位為MPa/nm)的速度分布曲線見圖3。

圖3 不同孔徑的密度(速度)分布曲線

由圖3 可見：不同孔徑內流體均為對稱分布，孔徑2 nm 狹縫內兩側各形成2 個吸附層，孔徑分別為5 nm 和8 nm 狹縫內兩側均形成四個吸附層，且各吸附層峰值密度逐漸減??；吸附相密度與狹縫孔徑大小有關，2 nm 孔徑中第一吸附層密度約為0.81 g/cm3，5 nm 孔徑和8 nm 孔徑中第一吸附層密度相近，約為0.92 g/cm3。由此可見，隨著孔徑寬度的增大，吸附層數量、第一吸附層密度均呈先增大后趨于穩定的趨勢。

值得注意的是，孔道邊界處流體受到壁面吸附作用以及z 方向的驅動力，而孔道中央處流體僅受到驅動力作用?？紫秲攘黧w分布及受力情況見圖4。

圖4 孔隙內流體分布及受力情況

由圖4 可見：孔徑分別為5 nm 和8 nm 狹縫內存在吸附相流體和體相流體，而2 nm 狹縫內充滿吸附相流體，不存在體相流體，這是由于在較小的孔隙中，無機質表面對烷烴分子的作用力強于烷烴分子內部的作用力，因此所有烷烴均吸附在無機質表面；當孔徑增大時，無機質表面對烷烴的引力逐漸減弱，因此在距離壁面較遠的孔道中央形成體相流體，呈游離態，密度為0.65 g/cm3，與NIST 提供的試驗數值一致。

根據圖3 中不同孔徑狹縫內流體的速度分布，其中散點為模擬結果，實線為該散點的拋物線擬合結果，孔徑分別為2，5，8 nm 狹縫內流體流速峰值分別為0.00029，0.00178，0.00599 ?/fs，結果表明在其他條件相同時，與小孔徑相比，大孔徑內流體流速遠大于小孔徑內的流體流速。狹縫內流體并非以同一速度流動，受到無機質壁面的影響，壁面附近流體流速低，孔道中央流體流速高，速度分布為拋物線形狀，可由二次函數擬合：

式中：v為流體流速，z為z方向x軸上距中心點的距離，a和c為擬合參數。

兩光滑平面內不可壓縮層狀流體的穩態速度剖面v(z)可通過Poiseuille 方程描述：

式中：n表示分子數，F表示驅動力，η表示流體黏度，W表示孔徑寬度。

根據連續流體力學理論假設，流體在固體壁面處速度為0，但烷烴在狹縫內流動時，固液界面存在滑移現象[8]，在納米級別體系中無法忽略同處在納米級別的滑移長度，因此引入滑移長度Ls建立更為通用的邊界條件，狹縫內具有滑移的Poiseuille 流動的速度剖面為：

烷烴黏度可表達為：

滑移長度定義為從無機質壁面到流體速度為0 處的外推距離[9]，在流體力學中，邊界條件模型包括正滑移模型、無滑移模型以及負滑移模型[10]，滑移長度模型見圖5。

圖5 滑移長度模型

圖5 中，vs為滑移速度，b 為滑移長度。b ＞0 時為正滑移模型，表示狹縫內所有流體均參與流動，且壁面處流體流速不為0；b=0 時為無滑移模型，表示壁面流體流速恰好為0；b＜0時為負滑移模型，表示壁面附近具有厚度為b 的黏滯層，該部分流體不參與流動。

不同孔徑的流體滑移長度及黏度見圖6。

圖6 不同孔徑的流體滑移長度及黏度

由圖6 可見：該驅動力條件下滑移長度均為負值，隨孔徑增大，負滑移長度先是急劇減小，后趨于平緩，這是由于小孔徑內兩側壁面距離很近，吸附在一側的流體仍會受到另一側壁面的吸引從而減小了單側壁面對流體的相互作用，隨著孔徑增大，受另一側壁面的影響減小，不參與流動的黏滯層厚度也因此增大。

針對流體黏度，分別對比狹縫內全部流體平均黏度與體相流體黏度。ηeff為平均黏度，ηcenter僅對中間體相流體擬合，為體相流體黏度。2 nm 孔徑內流體均為吸附相，不存在體相流體，因此ηeff和ηcenter基本一致?？讖椒謩e為5 nm 和8 nm 狹縫中存在體相流體，由于烷烴的剪切稀釋作用導致ηeff小于ηcenter?？讖綖?，5，8 nm 狹縫內流體的平均黏度分別為0.24，0.30，0.33 mPa·s，說明隨著孔徑的增大，流體的平均黏度有增大的趨勢，逐漸接近體相流體的黏度。

2.2 壓力梯度對流體運移影響

驅動力分別為0.002，0.0025，0.003 kcal/(mol·?)，孔徑為5 nm 狹縫內流體的速度分布曲線見圖7，不同壓力梯度下流體的滑移長度及黏度見圖8。

圖7 不同壓力梯度下流體的速度分布曲線

圖8 不同壓力梯度下流體的滑移長度及黏度

由圖7 和圖8 可見：不同驅替壓力梯度下孔道內流體的密度分布基本重合，因此受限空間內烷烴的賦存特征與驅動力無關。隨壓力梯度的增大，流體流動速度增大，滑移長度也增大，壓力梯度從0.002 kcal/(mol·?)增大到0.0025 kcal/(mol·?)時，滑移長度由-0.35 nm 增大到-0.13 nm，這是由于不流動的黏滯層中遠離壁面的部分逐漸被驅動，導致黏滯層厚度降低。壓力梯度從0.0025 kcal/(mol·?)增大到0.003 kcal/(mol·?)時，滑移長度由-0.13 nm 增大至0.004 nm，此時狹縫中不存在黏滯層，包括靠近壁面處的流體在內，狹縫內全部流體參與流動。在三種壓力梯度下，流體黏度幾乎不變，均在0.30 mPa·s 左右，可見狹縫內流體的平均黏度幾乎不受驅動力的影響；不同礦物孔隙中流體表觀黏度也有差異。

Zhang 等[11]研究了不同壓力梯度下方解石孔隙中辛烷黏度約為0.48 mPa·s，王森等[6]在石英孔隙中研究了不同壓力梯度下辛烷黏度約為0.29 mPa·s，與筆者的黏度結果十分接近，原因是混合無機孔隙中石英對烷烴的相互作用起到了較大的影響。

2.3 溫度對流體運移影響

不同溫度下第一吸附層峰值密度及體相流體密度見圖9。

圖9 第一吸附層峰值密度及體相流體密度

由圖9 可見：隨著溫度的升高，第一吸附層峰值密度由0.92 g/cm3變為0.81 g/cm3，有明顯降低趨勢，即升溫更有利于烷烴分子從無機質表面逃離，壁面對烷烴的吸附能力受到抑制，因此吸附相密度減小。不同溫度下流體的密度(速度)分布曲線見圖10，不同溫度下流體的滑移長度及黏度見圖11。

由圖10 和圖11 可見：升溫過程中，流體流速增大，平均黏度降低，流動有所增強；滑移長度由-0.35 nm 增大至0.8 nm，這意味著360 K 溫度下，狹縫兩側存在厚度為0.35 nm 的黏滯層；隨著溫度升高，滑移長度逐漸由負變為正，黏滯層厚度在逐漸降低，邊界處原本不流動的黏滯層開始慢慢被驅動，當滑移長度變為正值時，狹縫內全部流體均參與流動，因此熱采可成為提高頁巖油采收率的有效措施[12]。

2.4 含水率對流體運移影響

含水率分別為20%，50%，80%的油水平衡構型見圖12。

由圖12 可以直觀地看出當孔隙中存在油水兩相時，石英-伊利石混合壁面優先吸附水，含水較低時，在兩側壁面附近形成水膜，孔道中央為烷烴分子，形成水-油-水的夾層結構，隨著含水率的增大，壁面附近的水膜厚度增大，當含水率達80%時，混合無機納米孔內出現水橋，烷烴分子在孔道中央聚集成團簇狀。

不同的含水率條件下流體的密度及速度分布見圖13～15。

圖13 20%含水率時流體密度及速度分布

圖14 50%含水率時流體密度及速度分布

圖15 80%含水率時流體密度及速度分布

由圖13～15 可見：含水率為20%時，水在壁面處僅形成一個吸附層，吸附層峰值密度為1.101 g/cm3；當含水率增加到50%時，吸附層數量增加至4 個，第一吸附層峰值密度增大至1.26 g/cm3；含水率繼續增大至80%，吸附層數量和吸附層峰值密度趨于穩定，幾乎不再變化。這是由于含水率很低時，壁面吸附的水層并未達到飽和，隨著含水率的增加，吸附水層的數量、厚度也會隨之增加，當吸附達到飽和時，吸附水層的數量、厚度則趨于穩定。觀察烷烴的密度分布可以發現，含水率20%時，烷烴也存在微小的吸附效果，這是由于吸附水層厚度很小，不足以屏蔽壁面和烷烴間的相互作用，這種相互作用隨吸附水層的增大而減小，含水率50%時這種相互作用被水層完全屏蔽[13]。在含水率不斷增加的過程中，體相油區越來越小，油水兩相區域逐漸增大，當含水率達80%時，不再存在體相油區，此時邊界處為吸附水相，狹縫中其余空間均為油水混相區域。此時油水速度分布均呈拋物線形，Poiseuille模型擬合結果表明，油水速度分布呈明顯的不連續性，即烷烴和水之間存在液-液滑移現象。在第一吸附層附近水流速為0，這代表了混合無機納米孔表面的負滑移現象，存在一定厚度不參與流動的黏滯水層，滑移長度及油相黏度。不同含水率條件下烷烴滑移長度及油相黏度見表3。

表3 不同含水率時滑移長度及油相黏度

由表3 可見：不同含水率條件下烷烴黏度幾乎一致，可見油相黏度與油水比例無關。

3 結論

基于分子動力學方法研究了混合無機納米孔內流體的運移規律，通過非平衡分子模型對壓力驅動流動進行了模擬，討論了孔喉尺寸、壓力梯度、溫度以及含水率對流體賦存及運移的影響。對比不同條件下密度、速度、黏度及滑移長度等參數，得到以下結論。

1)孔徑表面會形成一定數量的吸附層，吸附層數量隨孔徑增大先增大后穩定，小孔徑狹縫內只存在吸附相流體，大孔徑狹縫內在邊界處形成吸附相流體，在孔道中央處形成游離態的體相流體。

2)相同壓力梯度下大孔徑內流體更容易被驅動， 增大壓力梯度以及升溫可以使邊界處的黏滯層也參與流動，同樣有利于納米孔中流體流動。

3)狹縫中存在油水兩相時，相對于烷烴，壁面優先吸附水分子，壁面附近存在水膜，含水率增大，水膜厚度先增大后趨于穩定，含水達80%時，烷烴成團簇狀聚集在孔道中央。

4)不同含水率時烷烴黏度幾乎一致，因此烷烴黏度與油水比例無關。