朱維耀，陳 震，尚新春

北京科技大學(xué)土木與資源工程學(xué)院，北京 100083

隨著常規(guī)油氣田的產(chǎn)能逐漸下降，非常規(guī)油氣藏的開發(fā)成為針對(duì)未來巨大能源需求供給的重要新興課題[1].在我國，非常規(guī)油氣藏的范圍包括頁巖油氣藏、致密油氣藏、煤層氣藏、油頁巖等多種資源類型，在國外還將稠油、油砂、天然氣水合物等也納入非常規(guī)油氣藏的范疇[2].這類油氣藏被廣泛開發(fā)主要?dú)w因于水平井鉆探和大規(guī)模水力壓裂（壓裂）技術(shù)的進(jìn)步[3].盡管已有不少成功的開采案例，但其產(chǎn)量的準(zhǔn)確預(yù)測仍具有很大的不確定性，這種不確定性源于多物理場耦合作用所導(dǎo)致的滲流場非線性特性加劇[4].因此，多場耦合問題是研究非常規(guī)油氣藏開采的一項(xiàng)重要科學(xué)問題.以往對(duì)多場耦合理論的運(yùn)用存在簡易化和適應(yīng)性缺陷，工程中缺少理論指導(dǎo)下的開采工藝與開發(fā)方法，制約著這類油氣藏的大規(guī)模高效開發(fā)，亟需對(duì)多場耦合滲流力學(xué)理論進(jìn)行深入認(rèn)識(shí).

油氣藏的多場耦合是指地下儲(chǔ)層中溫度場-滲流場-應(yīng)力場（Thermal-hydrological-mechanical,THM）的相互作用[5]，是巖層復(fù)雜條件中真實(shí)行為的抽象反映，是由固體力學(xué)、滲流力學(xué)、傳熱傳質(zhì)學(xué)等眾多工程學(xué)科相互交叉而形成的一門邊緣學(xué)科[6-9].非常規(guī)油氣藏開采過程中的多場耦合作用理論與工程問題的核心是不同物理場對(duì)非均質(zhì)、非線性滲流場的作用和影響[10]，對(duì)這種影響的實(shí)驗(yàn)驗(yàn)證、理論描述及數(shù)值模擬是研究中的重點(diǎn)[11]，具體包括應(yīng)力場變化對(duì)儲(chǔ)層滲流條件的影響、溫度場變化對(duì)滲流條件的影響、多場耦合作用對(duì)流體特性的影響，以及各物理場的理論模型之間如何進(jìn)行耦合和對(duì)應(yīng)的解耦、求解等問題，這些問題的解決都可為現(xiàn)場生產(chǎn)提供重要的理論、工程指導(dǎo)[12-15].

本文旨在了解非常規(guī)油氣藏開發(fā)過程中多場耦合問題的研究現(xiàn)狀與發(fā)展前景，提出了當(dāng)前多場耦合滲流力學(xué)理論發(fā)展所面臨的關(guān)鍵問題，針對(duì)數(shù)值模擬技術(shù)的現(xiàn)狀、進(jìn)展和瓶頸進(jìn)行了探討，并對(duì)多場耦合滲流力學(xué)理論及應(yīng)用的未來發(fā)展提出了展望.

1 非常規(guī)油氣開采多場耦合作用機(jī)制

非常規(guī)油氣藏的流體行為表現(xiàn)為多尺度多相流體（油、氣、水）在超低滲透、高度非均質(zhì)而且強(qiáng)應(yīng)力敏感的孔隙-裂縫介質(zhì)中流動(dòng).與常規(guī)油氣藏流動(dòng)相比，多場耦合作用更加顯著，造成高度非線性和復(fù)雜的流動(dòng)過程，并且伴隨劇烈的巖石-流體相互作用，以及納微米孔隙及微裂縫的變形[16-19].隨著巖石力學(xué)多物理場耦合試驗(yàn)以及各種現(xiàn)代化無損探測手段的發(fā)展，對(duì)非常規(guī)儲(chǔ)層在熱-流-固多場作用下的微觀演化過程和結(jié)構(gòu)變化的研究逐漸細(xì)化，可通過室內(nèi)實(shí)驗(yàn)量化地研究流固耦合作用、溫度的變化、有機(jī)質(zhì)特性的變化對(duì)流體流動(dòng)行為的影響規(guī)律等[20].

1.1 應(yīng)力場演化對(duì)滲流場的影響

應(yīng)力場演化對(duì)滲流場的影響通常采用巖石固有滲透率的變化來表征.在常規(guī)油氣藏中，應(yīng)力場的變化對(duì)巖石變形和流動(dòng)能力的影響一般較小，在實(shí)踐中大多會(huì)被忽略[20].但在非常規(guī)油藏中，應(yīng)力場的變化所引起的滲透率變化將顯著得多[21].通常認(rèn)為，開采過程中流體壓力的快速下降將導(dǎo)致有效應(yīng)力（總應(yīng)力與流體壓力之差）的上升，進(jìn)而使巖石進(jìn)一步被壓實(shí)并減小孔喉尺寸、減弱孔隙的連通性，因此降低滲透率.這一特性也被稱作滲透率應(yīng)力敏感性.研究者們通過將高溫高壓條件下的三軸試驗(yàn)系統(tǒng)與巖心流動(dòng)實(shí)驗(yàn)相結(jié)合，測試流固耦合作用下深部巖心的流動(dòng)能力[22]，并發(fā)現(xiàn)基質(zhì)巖心的表觀滲透率會(huì)隨著有效應(yīng)力的增大而降低，降低幅度可以達(dá)到30%～99%[23-27].一方面，黏土礦物的含量越高，應(yīng)力敏感性越顯著[28]；另一方面，非常規(guī)油氣藏中大量的烴類物質(zhì)吸附在有機(jī)質(zhì)的表面，一旦流體壓力降低，烴類物質(zhì)幾乎立即解吸[29]，導(dǎo)致吸附應(yīng)變[30]，宏觀表現(xiàn)為基質(zhì)骨架被壓縮、滲透率升高[31].不過，礦物含量和吸附對(duì)固體變形的影響相對(duì)較小，在大多數(shù)理論研究中都被忽略.

滲透率的應(yīng)力敏感性強(qiáng)度受儲(chǔ)層致密性的影響最大.筆者團(tuán)隊(duì)曾對(duì)巴西劈裂實(shí)驗(yàn)進(jìn)行改進(jìn)，利用單軸試驗(yàn)機(jī)緩慢施加載荷，可以在不破壞巖心整體構(gòu)型的情況下，得到人造微裂縫（圖1），從而對(duì)不同條件的非常規(guī)巖心進(jìn)行應(yīng)力敏感性測試[32].發(fā)現(xiàn)當(dāng)有效應(yīng)力從4 MPa 增至22 MPa 時(shí)，基質(zhì)巖心和裂縫巖心的滲透率下降幅度平均分別為95.5%和85%（圖2），即基質(zhì)的應(yīng)力敏感性相對(duì)更強(qiáng).為了在理論層面對(duì)基質(zhì)和裂縫進(jìn)行分別描述，多數(shù)研究者采用了多重介質(zhì)和離散裂縫模型[33].

圖1 裂縫性巖心的CT 掃描Fig. 1 CT scan image of a fractured core

圖2 基質(zhì)和裂縫巖心的應(yīng)力敏感性測試Fig. 2 Stress-sensitive experiment of matrix and fractured cores

影響滲透率應(yīng)力敏感性的主要因素還包括微觀孔道結(jié)構(gòu).對(duì)非常規(guī)儲(chǔ)層的壓裂改造將起到對(duì)微裂縫連通的作用，形成復(fù)雜的裂隙網(wǎng)絡(luò)，這使得應(yīng)力演化及應(yīng)力-應(yīng)變關(guān)系變得復(fù)雜，因?yàn)樵谟行?yīng)力的持續(xù)作用下，支撐劑可能被壓實(shí)，甚至可能發(fā)生壓碎，使裂縫閉合，導(dǎo)致不可逆的滲透率損失，巖石存在塑性變形行為[34-35].裂縫表面粗糙度越高，支撐能力越強(qiáng)，應(yīng)力敏感性越弱.不過，對(duì)地下真實(shí)裂縫的表面粗糙度進(jìn)行評(píng)估還是一項(xiàng)巨大的挑戰(zhàn)；而且，在數(shù)值模擬過程中應(yīng)用不同的應(yīng)力敏感數(shù)學(xué)模型將對(duì)數(shù)值計(jì)算結(jié)果產(chǎn)生明顯的影響，因此，如何從物理基本規(guī)律出發(fā)實(shí)現(xiàn)應(yīng)力敏感量化評(píng)價(jià)是亟待解決的問題.

值得注意的是，室內(nèi)的物理實(shí)驗(yàn)往往無法真實(shí)還原地下條件，這可能導(dǎo)致實(shí)驗(yàn)上的誤差.特別是非常規(guī)油氣藏的滲透率比常規(guī)儲(chǔ)層小3 到6 個(gè)數(shù)量級(jí)，這使得這類油氣藏的滲流實(shí)驗(yàn)非常耗時(shí)，也并非所有常規(guī)實(shí)驗(yàn)都適用于非常規(guī)儲(chǔ)層的測試[36-37].未來研究的重點(diǎn)應(yīng)在如何在實(shí)驗(yàn)室中精確模擬地下的真實(shí)應(yīng)力環(huán)境和流動(dòng)條件.

1.2 溫度場變化對(duì)滲流場的影響

非常規(guī)油氣藏儲(chǔ)層的溫度大都比常規(guī)儲(chǔ)層高，比如我國龍馬溪組頁巖氣藏的溫度可達(dá)150 ℃左右[38].在非常規(guī)油氣藏開采時(shí)必須進(jìn)行水力壓裂，這一過程可能會(huì)將多達(dá)兩萬至四萬噸水被注入地下[39]，與儲(chǔ)層發(fā)生能量交換.在常規(guī)油氣藏理論中很少考慮能量轉(zhuǎn)化效應(yīng)，但在非常規(guī)油氣藏開采過程中是必須考慮的.

首先，在非常規(guī)油氣藏中，巖石隨著溫度的變化而發(fā)生熱脹冷縮，這一效應(yīng)對(duì)流動(dòng)空間結(jié)構(gòu)的改變將是很顯著的[40]，固體骨架熱膨脹導(dǎo)致孔隙和裂縫系統(tǒng)空間的減少，從而導(dǎo)致滲透率的降低[41-42].其次，有機(jī)物是非常規(guī)油氣藏的重要組成部分，隨著熱成熟度的增加，有機(jī)質(zhì)發(fā)生熱解，析出原油、天然氣，增加巖石孔隙的體積以及流體壓力[43-44].還有研究表明，溫度上升所造成的烴類物質(zhì)解吸量增大也會(huì)導(dǎo)致巖石形變[45-47]，然而，烴類吸附為放熱過程、解吸為吸熱過程，因此發(fā)生解吸的位置儲(chǔ)層溫度將下降，反過來降低解吸量.研究人員必須通過分子動(dòng)力學(xué)等微觀方法計(jì)算宏觀量化指標(biāo)，如吸附熱[48]、吸附等溫線的交叉點(diǎn)[49]等.對(duì)烴類吸附以及置換過程中的熱量轉(zhuǎn)換和傳輸以及對(duì)物理場的影響機(jī)制還有待進(jìn)一步的量化研究.

此外，還必須考慮溫度的變化對(duì)流體黏度、密度的改變，特別是重質(zhì)油，其黏度及密度隨生產(chǎn)環(huán)境溫度變化的范圍可能很大[50].多種復(fù)雜的機(jī)理共同作用，導(dǎo)致目前對(duì)發(fā)生熱解后的巖石孔隙結(jié)構(gòu)的表征是非常規(guī)油氣藏開發(fā)的主要挑戰(zhàn)之一[51-52].現(xiàn)有研究通過一種或幾種方法表征部分尺度的孔隙空間隨溫度場的變化規(guī)律，如通過微觀CT 掃描技術(shù)來研究巖心加熱過程中的微米級(jí)裂縫和黏土礦物的演變[53-55]，通過掃描電鏡來觀察巖石孔隙的形態(tài)變化并劃分分布、類型[3,56]，以及通過氣體吸附實(shí)驗(yàn)對(duì)巖心樣品的孔隙結(jié)構(gòu)進(jìn)行分析[57].然而，因?yàn)槿狈τ行Ф煽康姆椒▉砭C合量化地比較使用不同方法獲得的多尺度數(shù)據(jù)，目前很少有關(guān)于溫度與多尺度孔隙之間綜合比較的研究.應(yīng)當(dāng)將這一問題作為下一步的研究方向，以期對(duì)溫度場的影響有更加全面的認(rèn)識(shí).

2 非常規(guī)油氣藏多場耦合作用理論模型

如上節(jié)所述，非常規(guī)油氣藏的儲(chǔ)層條件差，滲流場與應(yīng)力場、溫度場耦合作用，流體的流動(dòng)更為復(fù)雜，這就需要一套完整、精確的多場耦合作用理論，針對(duì)各個(gè)物理場建立多組偏微分方程（PDE），以通過不同流動(dòng)區(qū)域中的動(dòng)態(tài)物性模型來定義油氣的開采過程，從而區(qū)分不同區(qū)域之間以及多種物理場之間的相互作用.

2.1 應(yīng)力場的控制方程

滿足非常規(guī)油氣藏的應(yīng)力場控制方程可以通過擴(kuò)展已知的多孔介質(zhì)的本構(gòu)方程（應(yīng)變和應(yīng)力之間的關(guān)系）來獲得[58]，并由此得出固體質(zhì)點(diǎn)位移方程來描述介質(zhì)的變形.儲(chǔ)層機(jī)械行為通常被大多數(shù)研究假設(shè)為各向同性及彈性形變，再從基于Biot 原理的孔隙彈性理論和雙重介質(zhì)理論出發(fā)，通過基質(zhì)和裂縫的混合系統(tǒng)來描述機(jī)械變形的過程.當(dāng)頁巖或致密油氣藏投入開發(fā)后，孔隙及裂縫壓力降低，原始底應(yīng)力場發(fā)生改變，有效應(yīng)力增加，不僅使巖石骨架受到壓縮，也使固體質(zhì)點(diǎn)發(fā)生位移.通常近井地帶的位移較大，但隨著時(shí)間的增加，這種位移逐漸向遠(yuǎn)井地帶推移，這種隨時(shí)間變化的推移用固體靜力學(xué)變形方程進(jìn)行描述.被廣泛采用的油氣藏儲(chǔ)層變形控制方程如式(1)所示[59-61]：

其中，εij是總的應(yīng)變張量；G是巖石的剪切模量，Pa；KV是巖石的體積模量，Pa；σij表示巖體的應(yīng)力張量，Pa，i和j分別取1～3；σkk為巖體在三軸方向的法向應(yīng)力之和，Pa，σkk=σ11+σ22+σ33；δij為克羅內(nèi)克（Kronecker）函數(shù)，當(dāng)i=j時(shí)，函數(shù)取1，否則等于0；αm和αf分別為基質(zhì)系統(tǒng)和縫網(wǎng)系統(tǒng)的Biot系數(shù)，根據(jù)巖心力學(xué)實(shí)驗(yàn)測得[62-63]；pm和pf分別為基質(zhì)系統(tǒng)和縫網(wǎng)系統(tǒng)的孔隙壓力，Pa；εd為吸附引起的體積應(yīng)變.

剪切模量G和體積模量KV可通過下式計(jì)算[64]：

其中，E是巖石的彈性模量，Pa；ν是巖石的泊松比.

非常規(guī)油氣藏所特有的吸附-解吸行為應(yīng)在應(yīng)力場方程中加以考慮.通常，吸附被假設(shè)為僅發(fā)生在非常規(guī)油氣藏的基質(zhì)中，εd表示為如下的Langmuir 形式[61]：

其中，εL是Langmuir 吸附應(yīng)變；pL是Langmuir 壓力常數(shù)，Pa.

通常，開采過程中的固體力學(xué)問題被考慮為靜力學(xué)問題，此時(shí)滿足動(dòng)量守恒方程：σij,j+Fi=0.其中σij,j表示應(yīng)力 σij在j方向上的一階偏導(dǎo)數(shù).將動(dòng)量守恒方程以及應(yīng)變-位移幾何關(guān)系εi=引入式(1)，即轉(zhuǎn)化為以位移為待求量的偏微分方程組：

其中，ui,j代表位移ui（單位m）在j方向上的一階偏導(dǎo)數(shù)；ui,jk代表位移ui（單位m）在j和k方向上的二階偏導(dǎo)數(shù)；Fi為體積力，在力學(xué)問題中通常是物質(zhì)的重力，Pa.

當(dāng)考慮滲流場的影響時(shí)，應(yīng)力場控制方程(4)壓力pm及pf均下降，烴類物質(zhì)的解吸量增加，固體應(yīng)力的方程不再平衡，于是固體質(zhì)點(diǎn)發(fā)生位移，位移量的大小與壓降程度有關(guān)，位移的方向指向壓降最嚴(yán)重的區(qū)域.

當(dāng)考慮溫度場的影響時(shí)，在方程(4)中還應(yīng)增加線熱彈性形變項(xiàng)：KVαTTi，以表征溫度變化導(dǎo)致的固體熱脹冷縮.其中 αT為線熱彈性應(yīng)變系數(shù)，K-1；Ti為溫度T（單位K）在i方向上的一階偏導(dǎo)數(shù).

不過，上述模型中僅考慮了固體具有彈性應(yīng)變，而地下的巖石往往在具有彈性應(yīng)變的同時(shí)還具有塑性應(yīng)變，但目前考慮固體塑性應(yīng)變的模型并不多見，因?yàn)榈叵聨r石的塑形變形是很難觀測的.此外，應(yīng)力場的變化還不僅由于生產(chǎn)所導(dǎo)致，針對(duì)非常規(guī)油氣藏的水力壓裂以及生產(chǎn)過程中的重復(fù)壓裂都會(huì)使原有應(yīng)力環(huán)境發(fā)生改變，如使最大主應(yīng)力方向發(fā)生變化[65-66]、改變固體力學(xué)參數(shù)、裂縫延伸擴(kuò)展的模式發(fā)生轉(zhuǎn)變[67]等.這些復(fù)雜的因素都應(yīng)在未來的多場耦合研究中詳加考察并細(xì)致計(jì)算.

2.2 滲流場的控制方程

非常規(guī)油氣藏的滲流具有非線性特性，其滲流控制方程基于質(zhì)量守恒原理.多孔介質(zhì)（孔隙、裂縫）中多相流體（原油、天然氣、水）的存儲(chǔ)與流動(dòng)都用一般的滲流控制方程來表征：

其中，t表示時(shí)間，s；ρ表示流體在t時(shí)刻的密度，kg·m-3；v表示流體在t時(shí)刻的流速，m·s-1；q為質(zhì)量源匯項(xiàng)，kg·s-1；m表示單位多孔介質(zhì)體積中的流體質(zhì)量，kg，在非常規(guī)油氣藏中，此項(xiàng)包括游離相和吸附相[68]：

其中，ρ和 ρsc分別為烴類物質(zhì)在儲(chǔ)層和地面條件下的密度，kg·m-3；φ為多孔介質(zhì)巖體的孔隙度；ρr為多孔介質(zhì)巖體的密度，kg·m-3；p為流體壓力，MPa；VL為Langmuir 體積常數(shù)，m3·kg-1.

非常規(guī)油氣藏的非線性滲流特征體現(xiàn)在運(yùn)動(dòng)方程上.對(duì)于常規(guī)油氣藏，流體的速度v按達(dá)西定律的形式，表現(xiàn)為與有效滲透率K成正比的關(guān)系.然而，非常規(guī)油氣藏的流動(dòng)孔徑為納米級(jí)，流體流動(dòng)受納微米效應(yīng)的影響顯著，實(shí)測流量與線性達(dá)西定律（流速v與多孔介質(zhì)滲透率K成正比）[69]的理論值存在偏差，因此通常采用表觀滲透率Kapp（m2）來表征有效的滲透率[70].非常規(guī)油氣藏與常規(guī)儲(chǔ)層中的流動(dòng)條件和輸運(yùn)機(jī)理存在許多區(qū)別，滲流控制方程中必須整合處理流動(dòng)的主要機(jī)制，包括吸附、努森擴(kuò)散、滲流，以及溫度和壓力對(duì)流體物性的影響等[71-72].在一些合理假設(shè)下，裂縫及大孔隙中的流動(dòng)通過線性達(dá)西定律進(jìn)行描述，而納微米孔隙中的流動(dòng)則必須考慮各種流動(dòng)機(jī)制耦合影響[73-74]，在多場耦合條件下被采用得最多的流動(dòng)模型是Beskok 和Karniadakis[75]所建立的納微米通道傳輸方程（B-K 模型）：

其中，v為烴類物質(zhì)的流動(dòng)速度，m·s-1；μ為烴類物質(zhì)的黏度，Pa·s；Kr為相對(duì)滲透率；Kapp為烴類物質(zhì)流動(dòng)的表觀滲透率，m2；Keff為儲(chǔ)層變形時(shí)的固有滲透率，m2；α和b為實(shí)驗(yàn)室所測得的參數(shù)；Kn為多孔介質(zhì)流道的努森數(shù)，其表達(dá)式如式(8)所示[76]：

其中，λ為平均分子自由程，m；r為特征長度，在多孔介質(zhì)中為孔隙半徑，m；KB為玻爾茲曼常數(shù)，J·K-1；rc為分子碰撞有效直徑，當(dāng)考慮單一成分的分子時(shí)，等于分子直徑，m；φeff為儲(chǔ)層變形時(shí)的固有孔隙度；τ為多孔介質(zhì)的迂曲度.

從而，Kn的變化蛻化成為一個(gè)與φeff和Keff相關(guān)的量；其余的量可視為定值，或由實(shí)驗(yàn)數(shù)據(jù)測得其隨溫度、壓力的變化情況.方程(7)是對(duì)線性達(dá)西定律的非線性修正，一些研究報(bào)道了與方程(7)的形式不一致的流動(dòng)模型，這些流動(dòng)模型多采用滲流阻力代數(shù)和的形式，但也采用了表觀滲透率的表述形式，不過大多形式復(fù)雜、參數(shù)繁多，關(guān)鍵是大多僅有實(shí)驗(yàn)室尺度下的流動(dòng)模擬案例，而不能應(yīng)用于儲(chǔ)層宏觀理論的問題.具體可參閱Akilu 等[77]及Taghavinejad 等[78]的論文.

當(dāng)考慮應(yīng)力場的影響時(shí)，儲(chǔ)層骨架質(zhì)點(diǎn)發(fā)生位移，整體發(fā)生應(yīng)變，方程(6)中的孔隙度 φ隨應(yīng)變的增加而降低，方程(7)中的滲透率Keff隨孔隙度的降低而減小，從而流體流動(dòng)能力減弱、流量減小.因此，儲(chǔ)層降壓開采的同時(shí)流動(dòng)條件會(huì)變差.有多種模型描述這種介質(zhì)滲流能力的變化，見本文2.4 節(jié).

當(dāng)考慮溫度場的影響時(shí)，流體的密度ρ、黏度μ以及吸附量都將隨溫度的變化而改變.溫度越高，流體的密度、黏度越低，吸附量越小，因此滲流能力隨之增強(qiáng)、流量增大.反之亦然.

目前國內(nèi)外非常規(guī)油氣藏實(shí)際案例開始實(shí)施CO2置換等提高采收率辦法，但截至目前，仍缺乏一種行之有效地對(duì)混合流體的流動(dòng)行為進(jìn)行描述并得到充分驗(yàn)證的輸運(yùn)方程式.而且，非常規(guī)儲(chǔ)層中高溫、高壓條件下的流體相態(tài)變化也可能造成流動(dòng)條件的改變，盡管這一影響可能很重要，但目前還沒有一種綜合的方法來計(jì)算多組分流體在相變條件下的流動(dòng).這是未來值得研究的問題.

2.3 溫度場的控制方程

溫度場的演化可依據(jù)油氣藏能量方程來獲得.盡管同一個(gè)體積空間中共存著巖石骨架與流體，而且熱力學(xué)的特性不同，但依據(jù)溫度場作用普遍的基本假設(shè)：每個(gè)時(shí)間步中的熱量變化在瞬時(shí)完成，可將不同物質(zhì)的熱力學(xué)參數(shù)（如導(dǎo)熱系數(shù)、熱容等）作為等效的物理量，構(gòu)建能量守恒方程[79].其典型形式如式(9)所示：

其中，λc為巖石骨架與流體的等效導(dǎo)熱系數(shù)，W·m-1·K-1；T為溫度，K；Qst為單位摩爾吸附流體的吸附熱，J·mol-1；ρdsc為可吸附烴類物質(zhì)的地面條件密度，kg·m-3；Md為可吸附烴類物質(zhì)的摩爾質(zhì)量，kg·mol-1；Vd為當(dāng)前條件下單位體積巖石的吸附量，m3·kg-1；ηc為骨架與流體的質(zhì)量等效熱容，J·kg-1·K-1，其等效計(jì)算方法為：

其中，cr、co、cw、cg分別為巖石、油、水、氣的熱容，J·kg-1·K-1；So、Sw、Sg分別為油、水、氣的飽和度；ρo、ρw、ρg分別為油、水、氣在地下條件下的密度，kg·m-3.

通常，烴類化合物的吸附/解吸造成的熱量轉(zhuǎn)換使非常規(guī)儲(chǔ)層的溫度場發(fā)生變化的幅度并不大，僅有近井地帶由于解吸量較多，才導(dǎo)致發(fā)生一定的溫度降低[80]，通常可以忽略.但是，現(xiàn)在采用新型壓裂和開采工藝的案例逐漸受到關(guān)注.如無水液氮壓裂[81]、超臨界CO2壓裂[82]、電磁加熱開采[83]、微波加熱開采[84]等.這些新工藝和新技術(shù)在實(shí)施過程中，都面臨著顯著的溫度場大幅變化，通過多場耦合理論來準(zhǔn)確評(píng)估和預(yù)測這些過程是非常有必要的.

當(dāng)考慮滲流場的影響時(shí)，方程(9)中需要增加對(duì)流傳熱項(xiàng)，當(dāng)考慮油、水、氣三相流動(dòng)時(shí)，對(duì)流傳熱項(xiàng)為：

其中，vo、vw、vg分別為油、水、氣的流速，m·s-1.

該項(xiàng)保證了流體能借助傳輸過程將熱量傳遞到相鄰網(wǎng)格之中.流體的對(duì)流傳熱方向指向射孔位置，射孔位置是烴類解吸吸熱程度最大的位置，因此對(duì)流傳熱能夠降低儲(chǔ)層的溫度降低速率.

2.4 滲流-應(yīng)力-溫度三場耦合模型

對(duì)多種物理場相互作用、影響的模型進(jìn)行理論建模并聯(lián)立求解即為全耦合模型.方程(4)、(5)、(9)彼此獨(dú)立，聯(lián)立之后就構(gòu)成了非常規(guī)油氣藏中熱-流-固（THM）耦合模型的基本非線性物理過程，也是求解這一復(fù)雜過程的應(yīng)力、壓力以及溫度等物理量變化規(guī)律的基本方程形式.單獨(dú)的方程組較為理想化，各物理場之間的方程是獨(dú)立的，而一套完整的非常規(guī)油氣儲(chǔ)層的模擬器必須考慮各個(gè)物理場之間的參數(shù)影響，需要準(zhǔn)確關(guān)聯(lián)不同流動(dòng)和傳輸機(jī)制，將多物理場耦合在一起.物理場彼此之間的密切聯(lián)系如圖3 所示，整個(gè)過程是動(dòng)態(tài)且依賴時(shí)間的.

圖3 非常規(guī)油氣的多物理場之間的相互作用Fig. 3 Interactions between multiphysical fields in unconventional oil and gas reservoirs

（1）全耦合對(duì)滲流場的影響.

全耦合對(duì)滲流場的影響主要包括對(duì)滲透率和流體物性的影響.固體的滲透率受巖石形變的影響.根據(jù)巖石力學(xué)建立有關(guān)表征Keff的本構(gòu)方程，在多場耦合的研究中，被廣泛采用的本構(gòu)方程是三次法則[85]：

其中，K0是儲(chǔ)層未發(fā)生變形時(shí)的初始滲透率，m2；φ0分別為儲(chǔ)層未發(fā)生變形時(shí)的孔隙度.從而，Keff的變化問題又蛻化為φeff的變化問題.基于孔隙度變化量與固體應(yīng)變量的幾何關(guān)系，一些描述孔隙度-應(yīng)變的本構(gòu)模型已被提出并被用于應(yīng)力敏感巖石變化特性的實(shí)驗(yàn)驗(yàn)證[86-89].本團(tuán)隊(duì)博士亓倩對(duì)巖石骨架和孔隙的彈性模量進(jìn)行了研究，發(fā)現(xiàn)在前者遠(yuǎn)大于后者的情況下，基質(zhì)應(yīng)變遠(yuǎn)小于孔隙體應(yīng)變，可以僅考慮孔隙體變形，并建立了本構(gòu)模型[90]：

其中，αp為有效應(yīng)力系數(shù).

式(11)和(12)描述的是基質(zhì)內(nèi)的孔隙度、滲透率的變化.不過非常規(guī)油氣藏開采的關(guān)鍵是裂縫的導(dǎo)流能力，裂縫的表面粗糙度使裂縫能夠?qū)崿F(xiàn)自支撐，這種自支撐可以改變裂縫各向異性的滲透率[91]，這使裂縫介質(zhì)與基質(zhì)的演化存在區(qū)別.而且，非常規(guī)儲(chǔ)層的裂縫有效流動(dòng)寬度通常小于實(shí)際的裂縫寬度[92]，為裂縫應(yīng)力敏感物理模型的推導(dǎo)帶來難度.為此，現(xiàn)有研究大多采用經(jīng)驗(yàn)公式，如通過在各種巖性上對(duì)裂縫的加載和卸載行為進(jìn)行試驗(yàn)，從中得出裂縫的滲透率對(duì)應(yīng)力演化響應(yīng)的非線性關(guān)系[93].但是，現(xiàn)在還沒有有效的監(jiān)控地下裂縫表面粗糙度的手段，這是未來應(yīng)予考慮的問題.

流體物性受壓力和溫度的雙重作用，對(duì)其進(jìn)行數(shù)學(xué)描述是一項(xiàng)挑戰(zhàn).至今還沒有足夠精確同時(shí)方便計(jì)算調(diào)用的光滑的數(shù)學(xué)模型，只能通過經(jīng)驗(yàn)公式進(jìn)行描述，并在一定壓力、溫度范圍內(nèi)保證結(jié)果的精確性[94].在筆者的所做研究中，收集并調(diào)用了流體高壓高溫物性的數(shù)據(jù)庫，避免了使用任何經(jīng)驗(yàn)公式[80].采用這種方法比經(jīng)驗(yàn)公式準(zhǔn)確得多.

（2）應(yīng)力場和溫度場的相互作用.

應(yīng)力場和溫度場存在相互作用.主要體現(xiàn)在溫度的變化使固體產(chǎn)生熱脹冷縮的熱應(yīng)變，而固體的形變對(duì)溫度場產(chǎn)生應(yīng)變能.應(yīng)力場的控制方程中應(yīng)增加熱應(yīng)變項(xiàng)：

其中，εQ為熱應(yīng)變，Pa.

同時(shí)，溫度場的控制方程中應(yīng)根據(jù)應(yīng)變能的計(jì)算方法增加附加熱源：

其中，εV為當(dāng)前時(shí)間步中的固體應(yīng)變.

由此，才能形成滿足圖3 所示的方程耦合關(guān)系.這就是非常規(guī)油氣藏的全耦合基本模式.全耦合的計(jì)算需要聯(lián)立求解三個(gè)物理場的所有物理量，盡管準(zhǔn)確性高，編程非常復(fù)雜繁瑣、耗費(fèi)硬件資源，而且時(shí)間步長必須控制在很小的范圍內(nèi)以保證收斂.為此，研究人員們也通過半耦合模式來作為全耦合的替代.

2.5 滲流-應(yīng)力二場耦合模型

絕大多數(shù)情況下，非常規(guī)儲(chǔ)層僅采用天然能量開采，而不采用熱采或化學(xué)開采工藝，此時(shí)溫度場的影響相對(duì)較小，特定情況下可以忽略不計(jì)[80]，沒有必要計(jì)算滲流-應(yīng)力-溫度三場耦合，計(jì)算滲流-應(yīng)力二場耦合模型即可.因此，耦合模型僅包括滲流場和應(yīng)力場兩種物理場的控制方程，應(yīng)力的變化是儲(chǔ)層滲流變化的絕對(duì)主導(dǎo)因素，也被稱為儲(chǔ)層的應(yīng)力敏感性.

根據(jù)彈性力學(xué)，應(yīng)力場的控制方程通過有限元法進(jìn)行求解.但是，滲流場的控制方程大都利用有限差分或有限體積法求解.由于有限元法的網(wǎng)格與求解滲流場的網(wǎng)格（結(jié)構(gòu)化網(wǎng)格或非結(jié)構(gòu)化多邊形網(wǎng)格）不一致，所以在每個(gè)時(shí)間步計(jì)算過程中必須更換網(wǎng)格，所有數(shù)值在更換網(wǎng)格時(shí)均需要插值.為了避免頻繁更換網(wǎng)格所增加的不穩(wěn)定性，可選擇兩種解決方案.其一是滲流場也通過有限元法來求解，如使用商業(yè)軟件COMSOL；其二是簡化流固耦合過程，為此研究者們開發(fā)過許多經(jīng)驗(yàn)公式來描述滲透率隨壓力變化的關(guān)系，如最普遍的冪指數(shù)模型[61]：

其中，K0為初始滲透率，m2；c為基質(zhì)或裂縫的壓縮系數(shù)，Pa-1；p0為初始?jí)毫?，Pa；p為儲(chǔ)層的實(shí)時(shí)壓力，Pa.

由于非常規(guī)儲(chǔ)層裂縫系統(tǒng)的復(fù)雜性，有時(shí)經(jīng)驗(yàn)公式并不能完全擬合裂縫的滲透率或?qū)Я飨禂?shù)變化，但較為完整的實(shí)驗(yàn)數(shù)據(jù)可以作為數(shù)值計(jì)算的補(bǔ)充，模擬過程中對(duì)實(shí)驗(yàn)參數(shù)的直接調(diào)用往往是更加準(zhǔn)確的方法.在未來必須面臨的問題是對(duì)地下真實(shí)裂縫情況和行為的準(zhǔn)確認(rèn)識(shí)，特別是裂縫的自支撐作用對(duì)導(dǎo)流能力變化的影響，這可能與支撐劑的破碎性能及其在裂縫中的填充規(guī)律有關(guān)，但在目前的理論計(jì)算中這仍是一個(gè)巨大的挑戰(zhàn).

3 非常規(guī)油氣藏多場耦合作用仿真模擬

對(duì)多場耦合過程進(jìn)行模擬非常具有挑戰(zhàn)性：首先，開發(fā)出一個(gè)功能齊全的模擬器來模擬所有復(fù)雜物理過程是很難的；第二，找到一個(gè)能滿足理論上合理的解耦方法的本構(gòu)關(guān)系需要跨學(xué)科的專業(yè)知識(shí)，而且目前沒有通用的解決方案；第三，多場耦合仿真的模擬驗(yàn)證受到可以涵蓋所有物理場演化過程的嚴(yán)格限制.因而，問題的非線性很強(qiáng)，除了少數(shù)特殊情況可求出解析解外，一般需要用數(shù)值方法進(jìn)行求解.這就要求模擬方法滿足與實(shí)際情況吻合的準(zhǔn)確性、針對(duì)解耦過程的適應(yīng)性和可接受計(jì)算速度.

3.1 模擬仿真方法

流固耦合模擬仿真始于20 世紀(jì)80 年代，最初在單相流動(dòng)中得到實(shí)現(xiàn)[95].由于單相流的方程中少了相對(duì)滲透率、毛管壓力等一系列非線性因素，大多數(shù)仿真能夠以完全耦合的方式同時(shí)求解多物理場控制方程[96-97].多相流研究的重要性逐漸凸顯后，人們開始尋求不同的解耦方法，以提高計(jì)算效率和精度.這些耦合計(jì)算方法可分為單向耦合、迭代耦合和完全耦合三種.

（1）單向耦合（半耦合）.

模擬器在每個(gè)時(shí)間步中，先執(zhí)行滲流模擬計(jì)算，再將壓力和飽和度解傳遞給應(yīng)力場控制方程進(jìn)行應(yīng)力或位移的計(jì)算.這種方法也稱為松散耦合，因?yàn)閼?yīng)力場的計(jì)算只接收流動(dòng)解，而應(yīng)力或位移解并不反饋給滲流場的計(jì)算.考慮溫度場的情況也類似，即分別求解滲流場和應(yīng)力場、溫度場的三組方程[98-99]，也叫做半耦合.由于無需聯(lián)立求解所有物理場的參量，而是分步、依次地求解每個(gè)時(shí)間步內(nèi)的物理量，因此不是嚴(yán)格意義上的耦合，但是通過大量計(jì)算顯示，在時(shí)間步長得到有效控制的條件下，半耦合模型的結(jié)果也可以達(dá)到全耦合模型的精度[100]，因此也被研究者們采用.

（2）迭代耦合.

每個(gè)時(shí)間步中，分別依次求解各個(gè)物理場控制方程中的主要待求量，也被稱為順序耦合.通常先求解滲流場控制方程，然后將壓力和飽和度解傳給應(yīng)力場方程[101-102].應(yīng)力場方程的解再傳給到溫度場以及滲流場的控制方程，直到所有方程在同一個(gè)時(shí)間步中收斂[103].這種方法在每個(gè)時(shí)間步中都是雙向耦合的.這種方法實(shí)施起來較為簡便，占用內(nèi)存較小，運(yùn)算效率較高，穩(wěn)定性好，現(xiàn)今科研界被廣泛應(yīng)用的COMSOL 軟件即采用順序耦合.

（3）完全耦合.

將多個(gè)物理場的非線性偏微分方程在每個(gè)時(shí)間步中同時(shí)求解，一次性解得所有待求量.大多數(shù)完全耦合的仿真模擬僅限于單相流動(dòng)[104].這種方法的計(jì)算精度最高，不過處理和計(jì)算較難，而且對(duì)內(nèi)存和運(yùn)算能力的要求都很高，任何微小的初始化設(shè)置失誤都很容易導(dǎo)致計(jì)算不收斂.

除了上述常見的多場耦合研究外，還有一些研究者開發(fā)了適用于一定條件的半解析解耦方法.筆者基于非均勻全流場流固耦合參數(shù)的宏觀數(shù)學(xué)表征方法，離散單元上引入近似Boltzmann 變換，建立一種半數(shù)值-半解析求解多場耦合過程的方法，計(jì)算速度遠(yuǎn)高于數(shù)值求解，可在一階精度上解決非線性非穩(wěn)態(tài)非均勻的滲流問題.

3.2 考慮多場耦合作用的非常規(guī)油氣藏模擬

筆者在曾經(jīng)提出的三大區(qū)、五小區(qū)的裂縫網(wǎng)絡(luò)物理模型（圖4）上進(jìn)行多場耦合產(chǎn)能模擬研究，對(duì)于建立的非常規(guī)油氣藏的水平井體積壓裂非線性滲流偏微分方程及應(yīng)力場、溫度場的分區(qū)模型進(jìn)行構(gòu)建.目前的求解難點(diǎn)在于裂縫控制區(qū)縫網(wǎng)形態(tài)變化的數(shù)學(xué)描述，流固耦合計(jì)算過程需要對(duì)體積壓裂的復(fù)雜縫網(wǎng)結(jié)構(gòu)做理想化的處理，只能限制在彈性變形和較為簡單的彈塑性變形以內(nèi)，才能使用嚴(yán)格的數(shù)學(xué)公式對(duì)裂縫形態(tài)進(jìn)行描述.為此，筆者基于平板滲流理論、樹狀分形理論來表征體積壓裂復(fù)雜裂縫形態(tài)，分為樹狀、簇狀、羽狀、網(wǎng)狀四種基本形態(tài)，從而求取其等效滲透率以及應(yīng)力敏感變化[105].

圖4 “三大區(qū)，五小區(qū)”多尺度模型Fig. 4 Multiscale model of the “three main sectors,five subdivided sectors”

根據(jù)筆者所進(jìn)行的多相多場耦合計(jì)算產(chǎn)量變化曲線，溫度場變化主要受解吸附時(shí)溫度降低的影響，下降幅度較小，影響相對(duì)有限；而應(yīng)力場的變化對(duì)孔、滲的影響較大，流-固耦合作用對(duì)生產(chǎn)動(dòng)態(tài)有相對(duì)更顯著的影響.以非常規(guī)氣藏為例，考慮多場耦合作用時(shí)產(chǎn)氣量降低約13%（圖5）.因此地質(zhì)力學(xué)參數(shù)對(duì)產(chǎn)量的影響也較為顯著，如彈性模量越大，巖層抵抗變形的能力越強(qiáng)，表現(xiàn)出比較弱的多場耦合作用；即彈性模量越大，儲(chǔ)層剛度越大，流固耦合作用對(duì)流動(dòng)通道的影響越小，累產(chǎn)氣量、累產(chǎn)水量越大.并且，隨著地應(yīng)力的增加，巖石應(yīng)變?cè)鰪?qiáng)，孔隙度減小，產(chǎn)氣、產(chǎn)水均下降.若不考慮壓裂中應(yīng)力差對(duì)改造效果的影響，生產(chǎn)過程中氣井產(chǎn)能主要受控于最小水平主應(yīng)力大?。鹤钚∷街鲬?yīng)力越大，地層流體采出后，裂縫所承受的應(yīng)力傷害增加，累產(chǎn)氣降低.

圖5 多場耦合作用對(duì)多尺度流動(dòng)區(qū)域的影響Fig. 5 Multiphysical fields coupling effects on multiscale flow sectors

現(xiàn)代計(jì)算機(jī)運(yùn)算能力的飛速發(fā)展對(duì)通過數(shù)值方法高效和精確的求解提供了更大的可能性和前景[76].在闡述非常規(guī)油氣藏流動(dòng)和運(yùn)移機(jī)理的基礎(chǔ)上，相應(yīng)的數(shù)學(xué)關(guān)系可能被更加精細(xì)地描述并形成高效能的模擬器.盡管目前已在描述單個(gè)機(jī)制方面取得了進(jìn)展，但仍然缺乏能夠以綜合方式耦合大部分或所有機(jī)制的大型集成的油氣藏模擬器.另一方面，在普通油氣藏中考慮多物理場耦合問題會(huì)顯著增加模擬器所需的計(jì)算成本.為了提高非常規(guī)油氣藏模擬器的計(jì)算效率，需要研究可以進(jìn)行分布式計(jì)算或云計(jì)算的高性能計(jì)算方案，并專門用于大規(guī)模非常規(guī)油氣藏模擬.

4 結(jié)論與展望

對(duì)于非常規(guī)油氣藏的多物理場耦合滲流理論和數(shù)值模擬的研究，由單相、線性滲流理論向多相、非線性、跨尺度、多物理場耦合計(jì)算發(fā)展，特別是在近年來非常規(guī)資源愈加受到重視的時(shí)期，新發(fā)現(xiàn)的資源大多為含有納微米孔隙的非常規(guī)油氣藏，由于目前均采用體積壓裂后降壓開采，流場受到應(yīng)力場、溫度場作用明顯.建議對(duì)于室內(nèi)試驗(yàn)中模擬地下真實(shí)應(yīng)力環(huán)境、溫度條件、烴類吸附及置換的熱量變化等問題深入研究，深化對(duì)地下巖石的塑性應(yīng)變、重復(fù)壓裂過程中的應(yīng)力環(huán)境變化、混合烴類的輸運(yùn)模型、以及孔滲隨應(yīng)力和溫度發(fā)生變化的數(shù)學(xué)模型等問題的認(rèn)識(shí).依據(jù)所建立的適應(yīng)性多場耦合模型，配合研究分布式計(jì)算或云計(jì)算的高性能計(jì)算方案，用于大規(guī)模非常規(guī)油氣藏模擬.此外，筆者認(rèn)為針對(duì)多場耦合滲流理論的研究，在繼承前人研究結(jié)論認(rèn)識(shí)的基礎(chǔ)上，應(yīng)在理論深度、大型化計(jì)算等方面進(jìn)一步深化研究，使其推廣到更多的領(lǐng)域.具體展望如下：

（1）深海油氣資源開采：深水低溫環(huán)境、地層強(qiáng)度偏低，采用大規(guī)模叢式井開采，井下設(shè)備磨損快.導(dǎo)致水下生產(chǎn)系統(tǒng)難建設(shè)，開發(fā)方案難制定.再加上低溫條件下的烴類物質(zhì)易形成水合物，導(dǎo)致管道阻塞.亟需結(jié)合多場耦合滲流力學(xué)，形成深海油氣高效開發(fā)方法.

（2）水合物開采：水合物資源處在高壓、低溫、復(fù)雜相變的環(huán)境之中.不論是采用較為簡單的降壓開采、CO2置換，還是熱采、化學(xué)開采等方法，涉及水合物分解相變與氣-水在地層內(nèi)滲流傳熱等多個(gè)物理化學(xué)過程之間的耦合機(jī)理都十分復(fù)雜，亟需建立水合物開采多場耦合滲流理論.

（3）替代稠油火燒油層工藝技術(shù)：通常對(duì)于蒸汽熱采效果不好的稠油油藏，火燒油層是一種對(duì)稠油油藏的終極技術(shù)，會(huì)為儲(chǔ)層帶來許多不利影響，包括潛在的、不可逆轉(zhuǎn)的化學(xué)作用.該工藝進(jìn)行后，稠油儲(chǔ)層即廢置.亟需尋找可替代技術(shù)，高效開發(fā)熱采效果不良的稠油，包括近年來的化學(xué)劑降黏、CO2降黏技術(shù)等，這個(gè)過程中所涉及的多場耦合過程亟需加以研究.還需在熱物理化學(xué)方法上尋求變革性技術(shù)或新的突破技術(shù)，如核、控溫近火燒油層、微波技術(shù)的安全利用，以及結(jié)合立體井軌跡與流場多區(qū)域適配控制技術(shù)等，并建立更加前沿的適應(yīng)性多場耦合滲流力學(xué)理論.