氣體在雙重多孔介質(zhì)中運移的分?jǐn)?shù)階導(dǎo)數(shù)建模研究

胡 韜 吳本霞 余星吉 靳 勇

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

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

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

1 模型與方法介紹

1.1 分?jǐn)?shù)階對流—擴散模型

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

(1)

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

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

分?jǐn)?shù)階對流—擴散方程模型的表達式如下：

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

(2)

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

(3)

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

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

(ⅰ) Riemanm-Liouville定義[18]

(4)

和

(5)

(ⅱ) Caputo定義

(6)

(ⅲ) Grünwald-Letnikov定義

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

(7)

(8)

1.2 數(shù)值模擬方法

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

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

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

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

(9)

2.1 算例1

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

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

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

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

2.2 算例2

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

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

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

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

2.3 算例3

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

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

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

3 結(jié)果討論

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

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

圖4 利用時間分?jǐn)?shù)階對流—擴散方程(FCDE)模型和經(jīng)典對流—擴散方程(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)或流場中粒子運動的全域相關(guān)和歷史依賴性，已成為一種流行的建模工具。

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

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

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