胡 韜 吳本霞 余星吉 靳 勇

(南京市江南小化工集中整治工作現(xiàn)場(chǎng)指揮部，江蘇 南京 210038)

在我國(guó)經(jīng)濟(jì)快速發(fā)展和工業(yè)化進(jìn)程不斷前進(jìn)中，全球的能源如石油、煤和天然氣等占據(jù)重要地位，也是國(guó)民經(jīng)濟(jì)發(fā)展的重要基礎(chǔ)[1]。隨著油氣工業(yè)的快速發(fā)展以及油氣資源的不斷消耗，非常規(guī)油氣資源如致密油、致密砂巖氣、煤層氣和頁巖油氣等正逐步取代常規(guī)油氣資源[2]。目前，隨著自然資源日益減少，環(huán)境污染與能源開采問題是能源可持續(xù)發(fā)展的兩大關(guān)鍵問題[3]，越來越受到整個(gè)社會(huì)的高度重視。眾所周知，雙重多孔介質(zhì)系統(tǒng)如天然土壤、裂縫巖石、油氣儲(chǔ)藏、煤氣層等，其裂縫和基質(zhì)孔隙自身結(jié)構(gòu)具有不同的幾何特性。大多數(shù)油氣藏結(jié)構(gòu)也是由多孔介質(zhì)與樹狀分叉網(wǎng)絡(luò)鑲嵌而成的雙重多孔介質(zhì)結(jié)構(gòu)[4]。雙重多孔介質(zhì)中的流體流動(dòng)廣泛應(yīng)用于水文、地質(zhì)、石油工程、農(nóng)業(yè)工程、環(huán)境工程等自然界和工業(yè)生產(chǎn)中，如地下水源污染預(yù)測(cè)、二氧化碳地下埋存和油氣資源開采輸運(yùn)等，并越來越受到更多領(lǐng)域研究學(xué)者的高度關(guān)注和重視[5]。就當(dāng)前研究現(xiàn)狀而言，研究較多的是地下水、污染物及溶質(zhì)的運(yùn)移過程，可為水污染預(yù)測(cè)及治理提供理論依據(jù)[6]。而在油氣開采運(yùn)輸過程中，油氣儲(chǔ)藏存在分布范圍較廣、滲透性能較低和開采難度較大等特點(diǎn)，常用的注水驅(qū)替技術(shù)會(huì)損壞套管、發(fā)生水竄效應(yīng)且采收率較低。而氣驅(qū)法最早源于1958年美國(guó)進(jìn)行的注CO2混相驅(qū)項(xiàng)目；1996 年，美國(guó)又在圣胡安盆地 Burlington Allison試驗(yàn)區(qū)進(jìn)行了CO2驅(qū)替煤層氣現(xiàn)場(chǎng)試驗(yàn)，并達(dá)到了提高采收率的效果。1963年我國(guó)將該技術(shù)引入到大慶油田；2003 年我國(guó)在山西的沁水盆地南部進(jìn)行了CO2驅(qū)替煤層氣的單井試驗(yàn)。此外，高等[7]利用不同高溫高壓實(shí)驗(yàn)系統(tǒng)研究了CO2與原油相互作用的機(jī)理；黃等[8]建立了CO2驅(qū)替原油的多相滲流模型?？傊?，將CO2注入到油氣藏中，不僅可以實(shí)現(xiàn)CO2的地質(zhì)封存，減少大氣中CO2的排放量，也可以降低油的粘度且提高油和氣的產(chǎn)量，從而達(dá)到提高油氣采收率和改善環(huán)境的雙贏目的[9]。因此本文主要關(guān)注氣體在復(fù)雜油氣藏中的運(yùn)移過程，搞清氣體運(yùn)移規(guī)律，對(duì)油氣勘探、開發(fā)都有重要的指導(dǎo)意義。

雙重多孔介質(zhì)中的氣體運(yùn)移過程由通過大尺度裂縫的滲流流動(dòng)和基質(zhì)內(nèi)部的擴(kuò)散運(yùn)動(dòng)兩種不同的物理機(jī)制構(gòu)成[10]。Bear[11]第一個(gè)將對(duì)流—擴(kuò)散方程(Convection-Diffusion Equation, CDE)模型用于描述宏觀多孔介質(zhì)中溶質(zhì)的運(yùn)移過程。而在非均質(zhì)土壤或儲(chǔ)藏結(jié)構(gòu)中，由于流體自身特有的物理屬性和賦存壓力以及外界環(huán)境因素的干擾，使得流體流動(dòng)出現(xiàn)反常擴(kuò)散現(xiàn)象，即所測(cè)得的穿透曲線(Breakthrough Curves, BTCs)通常表現(xiàn)出早到達(dá)和明顯的長(zhǎng)拖尾現(xiàn)象[12]，CDE模型無法充分描述該現(xiàn)象。本文引入了分?jǐn)?shù)階導(dǎo)數(shù)模型的應(yīng)用，分?jǐn)?shù)階微分方程具有記憶性或非局部性，在工程、信息處理、物理、金融、流體等領(lǐng)域的應(yīng)用越來越廣泛，比整數(shù)階微分方程能更好地模擬自然物理過程和動(dòng)力系統(tǒng)過程[13]。2000年，Benson等[14]最早將分?jǐn)?shù)階導(dǎo)數(shù)模型應(yīng)用于地下水和溶質(zhì)運(yùn)移中，模擬Cape Cod實(shí)驗(yàn)場(chǎng)地中溴離子的反常擴(kuò)散現(xiàn)象；Meerschaert等[15]通過研究發(fā)現(xiàn)，初始階段污染物在土壤中運(yùn)移的BTCs會(huì)出現(xiàn)提前穿透，而在運(yùn)移后期會(huì)出現(xiàn)長(zhǎng)拖尾現(xiàn)象。此外，Ali等[16]建立了用于描述致密介質(zhì)中氣體運(yùn)移規(guī)律的時(shí)間分?jǐn)?shù)階對(duì)流—擴(kuò)散模型。El Amin等[17]提出了空間分?jǐn)?shù)階導(dǎo)數(shù)模型用于探討氣體運(yùn)移過程?？傊壳把芯扛咝ч_采油氣資源及突破技術(shù)瓶頸是學(xué)者們急需解決的科學(xué)和工程問題，而將力學(xué)、數(shù)學(xué)與石油工程、地球科學(xué)等學(xué)科實(shí)現(xiàn)深度交叉融合，可以更加有效地推動(dòng)非常規(guī)油氣資源的勘探開發(fā)。

1 模型與方法介紹

1.1 分?jǐn)?shù)階對(duì)流—擴(kuò)散模型

從以往的研究中可知，對(duì)流—擴(kuò)散方程可以準(zhǔn)確地揭示流體在均勻多孔介質(zhì)中的運(yùn)移機(jī)理和時(shí)間、空間對(duì)運(yùn)移過程的影響，同時(shí)還可以描述流體的反應(yīng)擴(kuò)散、質(zhì)量分布和熱量運(yùn)輸?shù)任锢憩F(xiàn)象，且形式簡(jiǎn)單便于計(jì)算。通常用于描述均勻介質(zhì)中的粒子運(yùn)移問題的基本方程即經(jīng)典對(duì)流—擴(kuò)散方程，它的表達(dá)式如式1所示。

(1)

式1中C為濃度，v為平均流速，K為擴(kuò)散系數(shù)。

在自然界和實(shí)際工程應(yīng)用中，由于復(fù)雜多孔介質(zhì)(土壤、裂隙巖體)結(jié)構(gòu)和物理性質(zhì)的非均勻性導(dǎo)致了反常擴(kuò)散現(xiàn)象的出現(xiàn)，利用CDE模型描述氣體運(yùn)移問題時(shí)存在不足。因?yàn)榱黧w或粒子的運(yùn)移過程是隨著時(shí)間、空間發(fā)生變化的，并不是恒定不變的。大量的實(shí)驗(yàn)觀測(cè)和數(shù)值模擬結(jié)果也顯示氣體運(yùn)移的濃度穿透曲線具有時(shí)間方向上的拖尾現(xiàn)象，和空間上的非高斯分布特征，而分?jǐn)?shù)階導(dǎo)數(shù)算子具有歷史記憶性和空間依賴性，因此引入了分?jǐn)?shù)階對(duì)流—擴(kuò)散方程模型。本文主要開展了雙重多孔介質(zhì)中氣體運(yùn)移機(jī)理、運(yùn)移模型、簡(jiǎn)化求解、模型參數(shù)、模型應(yīng)用以及運(yùn)移規(guī)律等研究。

分?jǐn)?shù)階對(duì)流—擴(kuò)散方程模型的表達(dá)式如下：

(ⅰ) 時(shí)間分?jǐn)?shù)階對(duì)流—擴(kuò)散模型

(2)

(ⅱ) 時(shí)空分?jǐn)?shù)階對(duì)流—擴(kuò)散模型

(3)

這里α為時(shí)間分?jǐn)?shù)階階數(shù)，β為空間分?jǐn)?shù)階階數(shù)。

分?jǐn)?shù)階導(dǎo)數(shù)的定義由時(shí)間或空間的非局部卷積算子表示，刻畫粒子反常擴(kuò)散過程。常用的分?jǐn)?shù)階定義包括Riemanm-Liouville、Caputo和Gründwald-Letnikov。下面簡(jiǎn)要介紹幾種常用的分?jǐn)?shù)階定義。

(ⅰ) Riemanm-Liouville定義[18]

(4)

和

(5)

(ⅱ) Caputo定義

(6)

(ⅲ) Grünwald-Letnikov定義

——法國(guó)農(nóng)業(yè)部稱：預(yù)計(jì)2018年法國(guó)葡萄酒總產(chǎn)量為46.1億升，比2017年增長(zhǎng)27%，超出過去5年平均產(chǎn)量達(dá)6%。在波爾多，盡管部分酒莊受到冰雹和霜霉病的影響，但預(yù)計(jì)產(chǎn)量還將達(dá)到歷史正常水平。

(7)

(8)

1.2 數(shù)值模擬方法

通常分?jǐn)?shù)階微分方程的解析解含有特殊函數(shù)或者復(fù)雜的級(jí)數(shù)，難以獲得方程的解析解，因而數(shù)值求解方法引起了學(xué)者們的廣泛關(guān)注。多種求解對(duì)流—擴(kuò)散方程的數(shù)值方法，如有限差分法(Finite Difference Method)、有限體積法(Finite Volume Method)、有限元方法(Finite Element Method)和無網(wǎng)格方法(Meshless Method)等被利用[19-21]。

本文主要采用的是較為靈活、簡(jiǎn)便的數(shù)值計(jì)算方法，即有限差分法求解分?jǐn)?shù)階微分方程，且求解分?jǐn)?shù)階偏微分方程具有適用范圍較廣、精度較高和穩(wěn)定性較好等特點(diǎn)。Liu等[22]建立了一種顯式差分格式求解時(shí)間分?jǐn)?shù)階擴(kuò)散方程，并證明了方法的穩(wěn)定性和收斂性。Meerschaert和Tadjeran[23]利用有限差分法有效地獲得了分?jǐn)?shù)階微分方程的數(shù)值解。此外，Chen等[24]也較好地將有限差分法應(yīng)用于求解分?jǐn)?shù)階對(duì)流—擴(kuò)散方程中，且方法的穩(wěn)定性和收斂性被有效地證明。丁志清[25]利用有限差分法求解了空間分?jǐn)?shù)階導(dǎo)數(shù)方程，并進(jìn)行了收斂性分析。本文采用隱式有限差分法建立數(shù)值求解方案，該方法的穩(wěn)定性和收斂性已被證明，這里不再贅述，感興趣的讀者請(qǐng)參考文獻(xiàn)[26-28]。

2 數(shù)值模擬結(jié)果分析

近年來的研究發(fā)現(xiàn)，整數(shù)階導(dǎo)數(shù)的非線性和變系數(shù)模型能夠描述某些特殊情況下的反常擴(kuò)散過程，但依然存在適用面窄、參數(shù)物理意義不清晰和長(zhǎng)時(shí)間歷程模擬不準(zhǔn)確等缺點(diǎn)[18]。本文引入分?jǐn)?shù)階導(dǎo)數(shù)模型，并對(duì)其進(jìn)行數(shù)值求解。將時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)進(jìn)行離散，可以寫成如下式9的形式[29]。

(9)

2.1 算例1

假設(shè)給定入口邊界處一個(gè)連續(xù)點(diǎn)源Cin且持續(xù)向模擬樣品入口端通入氣體，并在出口邊界定義了對(duì)流邊界條件，氣體濃度曲線呈不斷上升趨勢(shì)直至逐漸達(dá)到穩(wěn)定。給出初始濃度值為0，且出口邊界條件設(shè)置為自由邊界條件。

通過數(shù)值模擬方法對(duì)氣體粒子在一定介質(zhì)區(qū)域內(nèi)的運(yùn)移過程進(jìn)行數(shù)值模擬，如圖1所示。

圖1 算例1：已知初始條件和邊界條件下，給定不同分?jǐn)?shù)階階數(shù)α求解時(shí)間分?jǐn)?shù)階對(duì)流—擴(kuò)散方程(FCDE)模型在一定時(shí)間內(nèi)模擬獲得的氣體粒子的幾種無量綱濃度穿透曲線。模型參數(shù)如下：v=0.2，K=0.05，T=240和L=120。

在非均勻介質(zhì)結(jié)構(gòu)中，分?jǐn)?shù)階導(dǎo)數(shù)模型是模擬粒子反常擴(kuò)散的有效工具。圖1顯示，時(shí)間分?jǐn)?shù)階對(duì)流—擴(kuò)散模型模擬結(jié)果中的分?jǐn)?shù)階階數(shù)α越小，粒子運(yùn)移過程越緩慢且濃度達(dá)到最大值所需的時(shí)間越長(zhǎng)。結(jié)合模擬結(jié)果推斷，介質(zhì)結(jié)構(gòu)是影響濃度穿透曲線變化的主導(dǎo)因素，分?jǐn)?shù)階模型中的階數(shù)有效反映了介質(zhì)結(jié)構(gòu)對(duì)粒子運(yùn)移的作用。由于油氣儲(chǔ)藏基質(zhì)具有較強(qiáng)的非均質(zhì)性結(jié)構(gòu)，雙重介質(zhì)結(jié)構(gòu)中存在死孔隙，或部分運(yùn)移通道被阻塞、填充，這些都阻礙了氣體粒子的運(yùn)移過程。從圖1可以看出，隨著分?jǐn)?shù)階模型中的時(shí)間階數(shù)α變小，則顯示氣體運(yùn)移過程越慢。而α=1時(shí)，經(jīng)典對(duì)流—擴(kuò)散模型在達(dá)到最大濃度之前以最快的增長(zhǎng)趨勢(shì)變化，直至濃度達(dá)到峰值。

2.2 算例2

假定在入口邊界處給定一初始值為C0，而設(shè)置入口邊界條件為0，同時(shí)設(shè)置出口邊界為自由邊界條件。

分?jǐn)?shù)階導(dǎo)數(shù)模型的數(shù)值模擬結(jié)果展示如下。

圖2顯示在時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)模型中，通過數(shù)值方法求得幾種不同時(shí)間階數(shù)α下的濃度穿透曲線。在分?jǐn)?shù)階對(duì)流—擴(kuò)散模型中，當(dāng)0<α<1時(shí)，分?jǐn)?shù)階階數(shù)α對(duì)氣體運(yùn)移過程具有強(qiáng)有力的影響。在固定的控制位置1/5L處，隨著階數(shù)α的減小，濃度峰值越小，并且反常擴(kuò)散會(huì)出現(xiàn)緩慢的衰減。圖2中可以看出濃度穿透曲線的尾部斜率隨階數(shù)α的改變而改變，階數(shù)越小，則冪律拖尾現(xiàn)象越嚴(yán)重，拖尾就越長(zhǎng)，這是由于粒子運(yùn)移通道被阻塞且阻礙了氣體的運(yùn)動(dòng)。圖2(b)是半對(duì)數(shù)坐標(biāo)圖，可以觀察到不同時(shí)間階數(shù)α下，濃度穿透曲線的尾部斜率不同，說明了時(shí)間階數(shù)α是影響拖尾現(xiàn)象的主要參數(shù)。

圖2 算例2：已知初始條件和邊界條件，由時(shí)間分?jǐn)?shù)階對(duì)流—擴(kuò)散模型數(shù)值求解獲得的不同階數(shù)α下的濃度穿透曲線(a)。模型參數(shù)如下：v=0.2，K=0.01，T=200和L=100。(b)是半對(duì)數(shù)坐標(biāo)圖。

2.3 算例3

在邊界處給定一瞬時(shí)點(diǎn)源，利用時(shí)空分?jǐn)?shù)階對(duì)流—擴(kuò)散模型(3)模擬粒子運(yùn)移規(guī)律。初始值為C0，入口邊界值為0，右邊邊界條件為自由邊界。數(shù)值模擬結(jié)果如下。

從圖3可以得出，空間分?jǐn)?shù)階對(duì)流—擴(kuò)散模型能夠較好地模擬粒子的早到達(dá)現(xiàn)象。隨著空間階數(shù)的減小，早到達(dá)現(xiàn)象越明顯。結(jié)合實(shí)際情況分析，若粒子在雙重多孔介質(zhì)中運(yùn)移時(shí)，介質(zhì)中存在大尺度裂縫通道，促使粒子加速運(yùn)動(dòng)，因此出現(xiàn)了早期的超擴(kuò)散現(xiàn)象。氣體在雙重多孔介質(zhì)中的運(yùn)移過程與孔隙結(jié)構(gòu)和空間分布密切相關(guān)。簡(jiǎn)而言之，分?jǐn)?shù)階導(dǎo)數(shù)模型是描述具有記憶效應(yīng)的氣體運(yùn)移過程的有效工具，下面將通過實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證該模型。

圖3 展示的是時(shí)空分?jǐn)?shù)階對(duì)流—擴(kuò)散模型模擬的數(shù)值結(jié)果(a)，(b)圖是半對(duì)數(shù)下的坐標(biāo)圖。模型參數(shù)為：C0=20，v=0.1，K=0.05，T=200和L=20。

3 結(jié)果討論

基于文獻(xiàn)[30]中的一系列實(shí)驗(yàn)研究，探討不同條件下包括壓力、溫度、氣體流速和介質(zhì)不均勻性等，CO2驅(qū)替頁巖上CH4的動(dòng)力學(xué)理論機(jī)制和CO2驅(qū)替的有利條件。經(jīng)典對(duì)流—擴(kuò)散模型在模擬CO2在頁巖儲(chǔ)藏的裂隙和基質(zhì)系統(tǒng)中的傳質(zhì)過程時(shí)，存在明顯不足。因此借助分?jǐn)?shù)階導(dǎo)數(shù)模型研究氣體運(yùn)移規(guī)律，模擬氣體的濃度穿透曲線。

在非均勻頁巖儲(chǔ)藏中，CO2置換CH4的過程較為復(fù)雜，CO2的運(yùn)動(dòng)會(huì)出現(xiàn)反常擴(kuò)散過程，其本質(zhì)是時(shí)間上有記憶性和空間上有非局域性的過程，而整數(shù)階導(dǎo)數(shù)極限定義具有局域性不能準(zhǔn)確地描述這類反常擴(kuò)散過程。從圖4可得出，分?jǐn)?shù)階導(dǎo)數(shù)模型相比于傳統(tǒng)模型更具準(zhǔn)確性和實(shí)用性。分?jǐn)?shù)階對(duì)流—擴(kuò)散模型能夠較好地描述驅(qū)替頁巖過程中CO2在非均勻儲(chǔ)層中的次擴(kuò)散過程。分?jǐn)?shù)階導(dǎo)數(shù)模型因其具有記憶和遺傳、路徑依賴性，能夠較準(zhǔn)確地刻畫氣體粒子的反常擴(kuò)散過程。

圖4 利用時(shí)間分?jǐn)?shù)階對(duì)流—擴(kuò)散方程(FCDE)模型和經(jīng)典對(duì)流—擴(kuò)散方程(CDE)模型模擬CO2的濃度穿透曲線。模型參數(shù)分別為：L=50 cm，K=10-4 cm2/s，v=0.15 cm/s。

4 總結(jié)

分?jǐn)?shù)階導(dǎo)數(shù)是微分—積分形式表示的非局域性算子，能準(zhǔn)確地刻畫復(fù)雜介質(zhì)結(jié)構(gòu)或流場(chǎng)中粒子運(yùn)動(dòng)的全域相關(guān)和歷史依賴性，已成為一種流行的建模工具。

分?jǐn)?shù)階對(duì)流—擴(kuò)散方程模型在描述地下含水層、河流、土壤和裂隙巖體等復(fù)雜介質(zhì)中的流體擴(kuò)散或遷移過程時(shí)，具有參數(shù)少、物理意義清楚和描述準(zhǔn)確等優(yōu)點(diǎn)。

從數(shù)值模擬結(jié)果得出，分?jǐn)?shù)階導(dǎo)數(shù)模型能夠較準(zhǔn)確地刻畫粒子的反常擴(kuò)散現(xiàn)象，通過改變分?jǐn)?shù)階模型中的時(shí)間或空間階數(shù)，則可反映不同程度的次擴(kuò)散或超擴(kuò)散現(xiàn)象。

從實(shí)驗(yàn)擬合結(jié)果得出，整數(shù)階導(dǎo)數(shù)的局部極限定義不適合描述地質(zhì)儲(chǔ)存中氣體粒子運(yùn)動(dòng)過程具有的路徑和時(shí)間依賴特征。而分?jǐn)?shù)階導(dǎo)數(shù)模型能夠較好地?cái)M合實(shí)驗(yàn)數(shù)據(jù)，準(zhǔn)確刻畫了氣體的運(yùn)移過程。