吳 兵 陳 希 田 坤 喬向陽(yáng) 王永科 辛翠平

（1. 陜西延長(zhǎng)石油（集團(tuán)）有限責(zé)任公司研究院，陜西 西安 710075；2. 中國(guó)石油新疆油田公司工程技術(shù)研究院，新疆 克拉瑪依 834000）

0 引 言

了解致密多孔介質(zhì)中的滲流特征對(duì)非常規(guī)油藏的開(kāi)發(fā)具有重要意義。致密多孔介質(zhì)有著復(fù)雜的微觀結(jié)構(gòu)，孔隙表面與流體分子之間的相互作用導(dǎo)致流體滑動(dòng)和流體物理性質(zhì)不均勻分布。在潤(rùn)濕性的致密多孔介質(zhì)中，相互作用會(huì)引起負(fù)滑移現(xiàn)象，用滑移長(zhǎng)度表示?；崎L(zhǎng)度取決于邊界區(qū)域的厚度以及孔隙中黏度和密度的分布。以上機(jī)制導(dǎo)致致密多孔介質(zhì)中流體的滲流特征與經(jīng)典達(dá)西定律的偏離，并顯示出非線性滲流特征［1‐10］。

關(guān)于非線性滲流特性目前有3 類(lèi)模型：分段模型、多參數(shù)模型和分形模型。分段模型是最常見(jiàn)的模型，它將流動(dòng)階段分為3 個(gè)部分，非線性流動(dòng)階段由一條曲線描述，在試井和數(shù)值模擬中被廣泛應(yīng)用，但壓力梯度較小時(shí)模型不適用。多參數(shù)模型是從實(shí)驗(yàn)擬合結(jié)果獲得的不同系數(shù)，以描述非線性滲流階段［11‐14］；鄧英爾等［15］提出了壓力梯度的隱式滲流速度表達(dá)式，但是模型的某些參數(shù)沒(méi)有明確的物理意義；時(shí)宇等［16］提出了考慮邊界層效應(yīng)的兩參數(shù)模型，但沒(méi)有考慮實(shí)際的曲折流動(dòng)路徑。通過(guò)引入孔隙分形維數(shù)和曲折分形維數(shù)來(lái)考慮孔隙半徑分布和流動(dòng)路徑的曲折性，分形理論在描述致密多孔介質(zhì)孔隙結(jié)構(gòu)的復(fù)雜性有一定的優(yōu)勢(shì)。目前已經(jīng)建立了一些分形滲流模型，J.C.Cai［17］基于毛細(xì)管束模型提出了塑性流體的分形滲流模型，但是沒(méi)有考慮邊界層的影響，S.Huang 等［18］的模型引入了邊界層效應(yīng)，但沒(méi)有考慮滑移長(zhǎng)度和流體物理性質(zhì)的不均勻分布。

此前的研究提出的非線性滲流模型，能夠描述不同條件下的滲流特征，但沒(méi)有綜合考慮邊界層效應(yīng)和微觀效應(yīng)的影響。基于毛細(xì)管束模型，本文提出了致密多孔介質(zhì)中流體傳輸?shù)姆中文Ｐ?，可以?zhǔn)確表征致密多孔介質(zhì)的有效滲透率，為致密油藏的開(kāi)發(fā)和數(shù)值模擬提供參考。

1 模型建立

1.1 流體物理性質(zhì)的不均勻分布

當(dāng)涉及到微納米級(jí)孔隙時(shí)，無(wú)滑動(dòng)邊界不再適用，流體在孔隙表面會(huì)發(fā)生負(fù)滑移現(xiàn)象，即流體分子能以一定的速度運(yùn)動(dòng)?？紫侗砻娴幕扑俣扰c邊界區(qū)域的厚度有關(guān)，還與相互作用引起的流體黏度的分布有關(guān)。當(dāng)流體是水時(shí)，若孔隙表面為親水性，孔隙表面施加的吸引力將增加流動(dòng)阻力并降低滑動(dòng)速度；若孔隙表面為疏水性，孔隙表面和流體的相互作用變成排斥力，將減小流動(dòng)阻力并增加滑動(dòng)速度。滑移速度用滑移長(zhǎng)度來(lái)表征，即

式中：vs——滑移速度，m s；ls——滑移長(zhǎng)度，μm；v——沿流動(dòng)方向的速度，m s；z——沿徑向方向的距離，m；r0——某一具體的半徑，μm。

流體物理特性的非均質(zhì)分布可以簡(jiǎn)化為體相區(qū)域和邊界區(qū)域（圖1 （a）），曲折特性（圖1（b））將通過(guò)分形縮放理論在1.3 節(jié)中討論?？紫侗砻婧土黧w相互作用僅影響邊界區(qū)域的物理性質(zhì)，不影響體相區(qū)域的物理性質(zhì)。

圖1 實(shí)際流動(dòng)行為的簡(jiǎn)化和曲折表征Fig.1 Simplification of actual flow behavior and tortuous characterization

邊界區(qū)域中流體黏度與接觸角的相關(guān)表達(dá)式為

式中：μi——邊界區(qū)域黏度，mPa?s；μb——體相區(qū)域黏度，mPa?s；θ——接觸角，（°）。

滑動(dòng)長(zhǎng)度與邊界區(qū)域厚度、體相區(qū)域黏度和邊界區(qū)域黏度之間的表達(dá)式為

式中δ——邊界區(qū)域的厚度，μm。

流體密度與接觸角之間的關(guān)系為

式中：ρi——體相區(qū)域密度，kg/m3；ρb——邊界區(qū)域密度，kg/m3。

邊界區(qū)域的厚度可以通過(guò)實(shí)驗(yàn)方法或分子模擬來(lái)確定。在耗散粒子動(dòng)力學(xué)（DPD）模擬和實(shí)驗(yàn)結(jié)果的基礎(chǔ)上，X.F.Tian 等［19］提出了流體黏度、孔隙半徑和壓力梯度的綜合表達(dá)式為

式中：r——孔隙半徑，μm；h1，h2，h3——恒定值，h1= 0.257 63，h2= ?0.261，h3= ?0.419；μ—— 黏度，mPa?s；?p——壓力梯度，MPa m。

S.Huang 等［18］將Tian 的表達(dá)式修改為

其中h1=a1a3a4，h2=a2，h3=a4

式中：a1——黏度的影響系數(shù)，(mPa·s)?1；a2——喉部半徑的影響系數(shù)，μm?1；a3——壓力梯度的影響系數(shù)，(MPa m)?1；a4——量綱為1 的壓力梯度影響指數(shù)。

1.2 微納米級(jí)孔隙中的流體流動(dòng)

邊界區(qū)域和體相區(qū)域同時(shí)存在2 個(gè)速度場(chǎng)，并且微納米孔隙內(nèi)液體的速度分布是連續(xù)的。在邊界區(qū)域，由于孔隙表面的滑移現(xiàn)象，流體速度公式為

式中：vi(r)——邊界區(qū)域的流體速度，m s；p——流體壓力，MPa；l——流體流動(dòng)方向的距離，μm。

考慮不均勻的黏度和邊界區(qū)域的厚度，可以通過(guò)公式計(jì)算體相區(qū)域的速度分布，即

式中vb(r)——體相區(qū)域內(nèi)的流體速度，m s。

計(jì)算邊界區(qū)域和體相區(qū)域內(nèi)的流體質(zhì)量流量，對(duì)式（7）和式（8）進(jìn)行積分。

邊界區(qū)域貢獻(xiàn)的質(zhì)量流量公式為

式中Mi——邊界區(qū)域貢獻(xiàn)的質(zhì)量流量，kg/s。

體相區(qū)域貢獻(xiàn)的質(zhì)量流量公式為

式中Mb——體相區(qū)域貢獻(xiàn)的質(zhì)量流量，kg/s。

計(jì)算總質(zhì)量流量對(duì)式（9）和式（10）求和，即

式中M——總質(zhì)量流量，kg/s。

1.3 分形理論

如圖2 所示，多孔介質(zhì)孔隙的分布符合分形規(guī)律，累計(jì)的孔徑分布表達(dá)式為

圖2 實(shí)際巖心樣品中邊界區(qū)域和體相區(qū)域的存儲(chǔ)特征Fig.2 Storage characteristics of interfacial region and bulk region in actual core samples

式中：N(L>r)——孔隙半徑大于r的孔隙的數(shù)量；rmax——最大孔隙半徑，μm；Df——孔隙分形維數(shù)（本次采用二維分形維數(shù)）。

Df計(jì)算公式為

式中：d——?dú)W幾里得維數(shù)，二維中d= 2；?——孔隙度；rmin——最小孔隙半徑，μm。

得到半徑大于rmin的孔隙的總數(shù)為

對(duì)于二維空間，孔隙模型的橫截面面積計(jì)算公式為

式中A——橫截面面積，m2。

如圖1 所示，致密多孔介質(zhì)中的流動(dòng)路徑被認(rèn)為是具有不同橫截面尺寸的曲折管。由于流路的曲折性質(zhì)，實(shí)際長(zhǎng)度Lt通常大于特征長(zhǎng)度L0（圖2），根據(jù)自相似分形定律，實(shí)際長(zhǎng)度Lt為

式中：Lt( )r——實(shí)際長(zhǎng)度，m；L0——特征長(zhǎng)度，m；DT——曲折分形維數(shù)，DT越高，流動(dòng)路徑越曲折。

式中：τave——平均迂曲度；rave——平均孔隙半徑，μm。

1.4 致密多孔介質(zhì)中滲透率的分形模型

邊界區(qū)域的厚度通常很小，可以忽略二次項(xiàng)，三次冪項(xiàng)和第四次冪項(xiàng)，將邊界區(qū)域貢獻(xiàn)的質(zhì)量流量簡(jiǎn)化為

結(jié)合式（15）、式（24），微納米級(jí)孔隙的流體總質(zhì)量流量計(jì)算公式為

考慮到致密油藏流動(dòng)路徑的曲折性質(zhì)，并經(jīng)過(guò)進(jìn)一步的公式推導(dǎo)，得

式中Ke——致密多孔介質(zhì)的有效滲透率，10?3μm2。

有效滲透率與壓力梯度的厚度、孔隙分形維數(shù)、曲折分形維數(shù)、流體黏度和密度的異質(zhì)分布有關(guān)。根據(jù)式（2）、式（4），它也直接受接觸角的影響，即

B.M.Yu 等［20］提出了不考慮多種機(jī)制的分形滲透率模型，而S.Huang 等［18］提出了考慮邊界層效應(yīng)的分形滲透率。如不考慮多種機(jī)理，則將模型簡(jiǎn)化為Yu 模型。如僅考慮邊界層效應(yīng)，可以將模型簡(jiǎn)化為Huang 模型，即

2 模型驗(yàn)證

楊正明等［21］在進(jìn)行特低滲透油藏非線性流滲流數(shù)值模擬研究時(shí)，通過(guò)恒速壓汞實(shí)驗(yàn)獲得了具有超低滲透率的巖心樣品的孔喉分布。該實(shí)驗(yàn)中，大多數(shù)巖心最大喉道半徑小于6 μm，流體的黏度為1.005 mPa?s，完整的實(shí)驗(yàn)程序和高精度的實(shí)驗(yàn)設(shè)備保證了實(shí)驗(yàn)數(shù)據(jù)的高精度，且該實(shí)驗(yàn)數(shù)據(jù)已用于驗(yàn)證S.Huang［18］的滲透率模型。

本次模型的驗(yàn)證也采用上述實(shí)驗(yàn)數(shù)據(jù)，用有效滲透率Ke與絕對(duì)滲透率K的比值對(duì)模型進(jìn)行驗(yàn)證，如果模型計(jì)算結(jié)果能夠成功地捕獲實(shí)驗(yàn)數(shù)據(jù)，則可以驗(yàn)證模型的可靠性。選用5 塊不同氣測(cè)滲透率值的巖心實(shí)驗(yàn)數(shù)據(jù)，模型計(jì)算應(yīng)用的參數(shù)見(jiàn)表1。模型計(jì)算預(yù)測(cè)結(jié)果（線）與楊正明實(shí)驗(yàn)數(shù)據(jù)（點(diǎn)）的比較如圖3 所示，該模型的預(yù)測(cè)結(jié)果與實(shí)驗(yàn)數(shù)據(jù)一致，模型的結(jié)果與實(shí)驗(yàn)數(shù)據(jù)具有相似的變化趨勢(shì)，有效滲透率和絕對(duì)滲透率的比值隨著壓力梯度而增加。盡管預(yù)測(cè)結(jié)果與實(shí)驗(yàn)數(shù)據(jù)顯示出良好的相關(guān)性，但提供更多模型計(jì)算需要的參數(shù)（如孔隙度、接觸角、曲折分形維數(shù)、孔隙分形維數(shù)等），可以進(jìn)一步提高模型計(jì)算的準(zhǔn)確性。

表1 模型計(jì)算中應(yīng)用的相關(guān)參數(shù)Table.1 Related parameters applied in the model’s calculations

圖3 實(shí)驗(yàn)數(shù)據(jù)（點(diǎn)）與模型的預(yù)測(cè)結(jié)果（線）對(duì)比Fig.3 Comparison between data collected from experi‐ments(plots)and the model’s prediction results(lines)

3 參數(shù)敏感性

3.1 潤(rùn)濕性和壓力梯度

圖4 為不同接觸角對(duì)有效滲透率的影響情況，不同接觸角下的計(jì)算結(jié)果是本次模型計(jì)算結(jié)果，Huang 模型和Yu 模型均未考慮到接觸角。因?yàn)镠uang 模型中邊界層效應(yīng)僅被視為有效孔隙半徑和質(zhì)量流量的減小，故本次模型計(jì)算的結(jié)果比Huang模型高。Yu 模型忽略了邊界層效應(yīng)和微觀效應(yīng)，故本次模型和Huang 模型的結(jié)果都比Yu 模型小。

圖4 不同接觸角對(duì)有效滲透率的影響Fig.4 Effects of contact angle on effective permeability

有效滲透率隨著壓力梯度的增加而增加，當(dāng)壓力梯度較大時(shí)趨于穩(wěn)定。這是由于壓力梯度的增加會(huì)導(dǎo)致邊界區(qū)域的厚度變小，滑移長(zhǎng)度增加；當(dāng)壓力梯度非常小時(shí)，流體不能克服由孔隙表面和流體分子之間的相互作用力，導(dǎo)致有效滲透率非常小。隨著壓力梯度的增加，邊界區(qū)域的厚度逐漸穩(wěn)定，有效滲透率的值也趨于恒定。

有效滲透率隨著接觸角的增大而增大，這是因?yàn)榻佑|角越大，邊界區(qū)域和體相區(qū)域中流體黏度差越小，滑移長(zhǎng)度越大，由孔隙表面引起的吸引力越小，有效滲透率值越大。故潤(rùn)濕性在致密多孔介質(zhì)中的液體傳輸能力中起著至關(guān)重要的作用，改善潤(rùn)濕性是提高致密油藏產(chǎn)量的有效措施。

3.2 流體黏度

如圖5 所示，流體黏度對(duì)有效滲透率有顯著影響，并且可以觀察到啟動(dòng)壓力梯度。有效滲透率隨著流體黏度的增加而降低，這是因?yàn)轲ざ仍酱?，邊界區(qū)域厚度越大，導(dǎo)致滲透率的降低。

圖5 流體黏度對(duì)有效滲透率的影響Fig.5 Effects of fluid viscosity on effective permeability

如圖6 所示，啟動(dòng)壓力梯度存在并且隨流體黏度而變化。一定的壓力梯度下，有效滲透率隨流體黏度的增加而減小。當(dāng)流體黏度達(dá)到一定程度時(shí)，壓力梯度較低時(shí)流體無(wú)法克服表面與液體之間的相互作用力。增加壓力梯度可以提高致密多孔介質(zhì)中的流體流動(dòng)能力。

圖6 壓力梯度對(duì)有效滲透率的影響Fig.6 Effects of pressure gradient on effective permeability

3.3 孔隙分形維數(shù)和曲折分形維數(shù)

從圖7 中可以看出，有效滲透率隨著孔隙分形維數(shù)的增大而增大。因?yàn)檩^大的分形維數(shù)意味著孔隙數(shù)增加，孔隙表面與流體之間相互作用的影響雖然更加明顯，但越來(lái)越多的孔隙仍然可以為流體提供傳輸路徑，所以當(dāng)孔隙分形維數(shù)增加時(shí)，壓力梯度的影響將變大。

圖7 孔隙分形維數(shù)對(duì)有效滲透率的影響Fig.7 Effects of pore fractal dimension on effective permeability

圖8 顯示有效滲透率隨曲折分形維數(shù)的增加而降低，這是因?yàn)榍鄯中尉S數(shù)越高，意味著曲折的流路越多，由于長(zhǎng)流路的阻力越大，有效滲透率就越低。當(dāng)曲折分形維數(shù)較大（DT>1.6）時(shí)，壓力梯度的影響會(huì)降低。

圖8 曲折分形維數(shù)對(duì)有效滲透率的影響Fig.8 Effects of tortuosity fractal dimension on effective permeability

4 結(jié) 論

（1）通過(guò)考慮邊界層效應(yīng)、負(fù)滑移現(xiàn)象和流體物理性質(zhì)不均勻分布等微觀效應(yīng)，基于毛細(xì)管束模型和分形理論，建立了新的滲透率模型來(lái)表征致密多孔介質(zhì)中流體的傳輸能力，并用高精度實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證了模型的可靠性。

（2）致密多孔介質(zhì)中流體有效滲透率隨著壓力梯度的增大而增大，隨著接觸角的增大而增大，隨著黏度的增大而減小；有效滲透率直接受到微觀孔隙結(jié)構(gòu)的影響，與孔隙分形維數(shù)呈正相關(guān)，與曲折分形維數(shù)呈負(fù)相關(guān)。