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

(1. 東北石油大學(xué)石油工程學(xué)院，黑龍江大慶 163318；2. 東北石油大學(xué)提高油氣采收率教育部重點(diǎn)實(shí)驗(yàn)室，黑龍江大慶 163318；3. 國(guó)家管網(wǎng)集團(tuán)油氣調(diào)控中心，北京 100013；4. 大慶油田第八采油廠，黑龍江大慶 163514；5. 大慶油田第四采油廠，黑龍江大慶 163511)

隨著社會(huì)經(jīng)濟(jì)發(fā)展，對(duì)油氣資源的需求日益增大，非常規(guī)油氣資源越來(lái)越受到重視，頁(yè)巖油也隨之成為世界各國(guó)油氣勘探的焦點(diǎn)。早在2012 年，美國(guó)便依靠大力開(kāi)發(fā)頁(yè)巖油實(shí)現(xiàn)了原油產(chǎn)量劇增，2019 年美國(guó)超過(guò)俄羅斯和沙特阿拉伯成為全球最大的原油生產(chǎn)國(guó)。由此可見(jiàn)，頁(yè)巖油的勘探開(kāi)發(fā)對(duì)世界能源格局有著不可忽視的影響[1]。

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

筆者以辛烷和水為研究對(duì)象，采用分子動(dòng)力學(xué)模擬方法開(kāi)展了石英和伊利石混合的無(wú)機(jī)質(zhì)孔隙內(nèi)流體賦存狀態(tài)和運(yùn)移規(guī)律的研究，通過(guò)分析孔隙內(nèi)流體的密度分布、速度分布、滑移長(zhǎng)度、黏度等參數(shù)，分別討論了孔徑大小、溫度、壓力梯度及含水率對(duì)納米孔隙內(nèi)流體賦存和運(yùn)移的影響。

1 模型構(gòu)建及模擬方案

1.1 模型構(gòu)建

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

表1 伊利石模型原子坐標(biāo)

表2 伊利石模型晶胞參數(shù)

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

圖1 伊利石-二氧化硅組合模型結(jié)構(gòu)

圖2 伊利石-二氧化硅混合頁(yè)巖孔隙內(nèi)流體分布初始結(jié)構(gòu)

1.2 模擬方案

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

2 結(jié)果分析

2.1 孔喉尺寸對(duì)流體運(yùn)移影響

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

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

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

值得注意的是，孔道邊界處流體受到壁面吸附作用以及z 方向的驅(qū)動(dòng)力，而孔道中央處流體僅受到驅(qū)動(dòng)力作用?？紫秲?nèi)流體分布及受力情況見(jiàn)圖4。

圖4 孔隙內(nèi)流體分布及受力情況

由圖4 可見(jiàn)：孔徑分別為5 nm 和8 nm 狹縫內(nèi)存在吸附相流體和體相流體，而2 nm 狹縫內(nèi)充滿吸附相流體，不存在體相流體，這是由于在較小的孔隙中，無(wú)機(jī)質(zhì)表面對(duì)烷烴分子的作用力強(qiáng)于烷烴分子內(nèi)部的作用力，因此所有烷烴均吸附在無(wú)機(jī)質(zhì)表面；當(dāng)孔徑增大時(shí)，無(wú)機(jī)質(zhì)表面對(duì)烷烴的引力逐漸減弱，因此在距離壁面較遠(yuǎn)的孔道中央形成體相流體，呈游離態(tài)，密度為0.65 g/cm3，與NIST 提供的試驗(yàn)數(shù)值一致。

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

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

兩光滑平面內(nèi)不可壓縮層狀流體的穩(wěn)態(tài)速度剖面v(z)可通過(guò)Poiseuille 方程描述：

式中：n表示分子數(shù)，F(xiàn)表示驅(qū)動(dòng)力，η表示流體黏度，W表示孔徑寬度。

根據(jù)連續(xù)流體力學(xué)理論假設(shè)，流體在固體壁面處速度為0，但烷烴在狹縫內(nèi)流動(dòng)時(shí)，固液界面存在滑移現(xiàn)象[8]，在納米級(jí)別體系中無(wú)法忽略同處在納米級(jí)別的滑移長(zhǎng)度，因此引入滑移長(zhǎng)度Ls建立更為通用的邊界條件，狹縫內(nèi)具有滑移的Poiseuille 流動(dòng)的速度剖面為：

烷烴黏度可表達(dá)為：

滑移長(zhǎng)度定義為從無(wú)機(jī)質(zhì)壁面到流體速度為0 處的外推距離[9]，在流體力學(xué)中，邊界條件模型包括正滑移模型、無(wú)滑移模型以及負(fù)滑移模型[10]，滑移長(zhǎng)度模型見(jiàn)圖5。

圖5 滑移長(zhǎng)度模型

圖5 中，vs為滑移速度，b 為滑移長(zhǎng)度。b ＞0 時(shí)為正滑移模型，表示狹縫內(nèi)所有流體均參與流動(dòng)，且壁面處流體流速不為0；b=0 時(shí)為無(wú)滑移模型，表示壁面流體流速恰好為0；b＜0時(shí)為負(fù)滑移模型，表示壁面附近具有厚度為b 的黏滯層，該部分流體不參與流動(dòng)。

不同孔徑的流體滑移長(zhǎng)度及黏度見(jiàn)圖6。

圖6 不同孔徑的流體滑移長(zhǎng)度及黏度

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

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

2.2 壓力梯度對(duì)流體運(yùn)移影響

驅(qū)動(dòng)力分別為0.002，0.0025，0.003 kcal/(mol·?)，孔徑為5 nm 狹縫內(nèi)流體的速度分布曲線見(jiàn)圖7，不同壓力梯度下流體的滑移長(zhǎng)度及黏度見(jiàn)圖8。

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

圖8 不同壓力梯度下流體的滑移長(zhǎng)度及黏度

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

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

2.3 溫度對(duì)流體運(yùn)移影響

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

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

由圖9 可見(jiàn)：隨著溫度的升高，第一吸附層峰值密度由0.92 g/cm3變?yōu)?.81 g/cm3，有明顯降低趨勢(shì)，即升溫更有利于烷烴分子從無(wú)機(jī)質(zhì)表面逃離，壁面對(duì)烷烴的吸附能力受到抑制，因此吸附相密度減小。不同溫度下流體的密度(速度)分布曲線見(jiàn)圖10，不同溫度下流體的滑移長(zhǎng)度及黏度見(jiàn)圖11。

由圖10 和圖11 可見(jiàn)：升溫過(guò)程中，流體流速增大，平均黏度降低，流動(dòng)有所增強(qiáng)；滑移長(zhǎng)度由-0.35 nm 增大至0.8 nm，這意味著360 K 溫度下，狹縫兩側(cè)存在厚度為0.35 nm 的黏滯層；隨著溫度升高，滑移長(zhǎng)度逐漸由負(fù)變?yōu)檎?，黏滯層厚度在逐漸降低，邊界處原本不流動(dòng)的黏滯層開(kāi)始慢慢被驅(qū)動(dòng)，當(dāng)滑移長(zhǎng)度變?yōu)檎禃r(shí)，狹縫內(nèi)全部流體均參與流動(dòng)，因此熱采可成為提高頁(yè)巖油采收率的有效措施[12]。

2.4 含水率對(duì)流體運(yùn)移影響

含水率分別為20%，50%，80%的油水平衡構(gòu)型見(jiàn)圖12。

由圖12 可以直觀地看出當(dāng)孔隙中存在油水兩相時(shí)，石英-伊利石混合壁面優(yōu)先吸附水，含水較低時(shí)，在兩側(cè)壁面附近形成水膜，孔道中央為烷烴分子，形成水-油-水的夾層結(jié)構(gòu)，隨著含水率的增大，壁面附近的水膜厚度增大，當(dāng)含水率達(dá)80%時(shí)，混合無(wú)機(jī)納米孔內(nèi)出現(xiàn)水橋，烷烴分子在孔道中央聚集成團(tuán)簇狀。

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

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

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

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

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

表3 不同含水率時(shí)滑移長(zhǎng)度及油相黏度

由表3 可見(jiàn)：不同含水率條件下烷烴黏度幾乎一致，可見(jiàn)油相黏度與油水比例無(wú)關(guān)。

3 結(jié)論

基于分子動(dòng)力學(xué)方法研究了混合無(wú)機(jī)納米孔內(nèi)流體的運(yùn)移規(guī)律，通過(guò)非平衡分子模型對(duì)壓力驅(qū)動(dòng)流動(dòng)進(jìn)行了模擬，討論了孔喉尺寸、壓力梯度、溫度以及含水率對(duì)流體賦存及運(yùn)移的影響。對(duì)比不同條件下密度、速度、黏度及滑移長(zhǎng)度等參數(shù)，得到以下結(jié)論。

1)孔徑表面會(huì)形成一定數(shù)量的吸附層，吸附層數(shù)量隨孔徑增大先增大后穩(wěn)定，小孔徑狹縫內(nèi)只存在吸附相流體，大孔徑狹縫內(nèi)在邊界處形成吸附相流體，在孔道中央處形成游離態(tài)的體相流體。

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

3)狹縫中存在油水兩相時(shí)，相對(duì)于烷烴，壁面優(yōu)先吸附水分子，壁面附近存在水膜，含水率增大，水膜厚度先增大后趨于穩(wěn)定，含水達(dá)80%時(shí)，烷烴成團(tuán)簇狀聚集在孔道中央。

4)不同含水率時(shí)烷烴黏度幾乎一致，因此烷烴黏度與油水比例無(wú)關(guān)。