朱維耀，陳 震，尚新春

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

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

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

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

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

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

1.1 應力場演化對滲流場的影響

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

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

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

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

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

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

1.2 溫度場變化對滲流場的影響

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

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

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

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

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

2.1 應力場的控制方程

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

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

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

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

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

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

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

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

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

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

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

2.2 滲流場的控制方程

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

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

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

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

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

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

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

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

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

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

2.3 溫度場的控制方程

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

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

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

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

當考慮滲流場的影響時，方程(9)中需要增加對流傳熱項，當考慮油、水、氣三相流動時，對流傳熱項為：

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

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

2.4 滲流-應力-溫度三場耦合模型

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

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

（1）全耦合對滲流場的影響.

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

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

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

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

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

（2）應力場和溫度場的相互作用.

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

其中，εQ為熱應變，Pa.

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

其中，εV為當前時間步中的固體應變.

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

2.5 滲流-應力二場耦合模型

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

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

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

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

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

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

3.1 模擬仿真方法

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

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

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

（2）迭代耦合.

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

（3）完全耦合.

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

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

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

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

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

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

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

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

4 結(jié)論與展望

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

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

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

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