曹靜，趙輝，喻高明，胡文明

(1.長(zhǎng)江大學(xué)石油工程學(xué)院，湖北武漢430100；2.中原油田普光氣田凈化廠，四川達(dá)州636100)

油藏?cái)?shù)值模擬是油藏開(kāi)發(fā)研究的重要手段之一.它能使工程師更好地了解油藏的各種物理性質(zhì)和流體在其中的流動(dòng)規(guī)律，以便作出正確評(píng)價(jià)，確定合理開(kāi)發(fā)方案和提高采收率.傳統(tǒng)油藏模擬器將運(yùn)動(dòng)方程、狀態(tài)方程和連續(xù)方程構(gòu)成的偏微分方程組離散化，轉(zhuǎn)化成一組非線性代數(shù)方程組，再采用迭代法求解.然而，由于實(shí)際油田建立的油藏模型可能包括幾十萬(wàn)或幾百萬(wàn)個(gè)網(wǎng)格，需要求解的方程數(shù)目非常大，所以，傳統(tǒng)油藏模擬器的數(shù)值解決方法非常耗時(shí).特別的，當(dāng)把油藏模擬器應(yīng)用于閉環(huán)油藏管理時(shí)[1?5]，其中的生產(chǎn)優(yōu)化和歷史擬合過(guò)程都需要反復(fù)運(yùn)行模擬器，導(dǎo)致它的計(jì)算成本大大提高.因此，在保證有足夠精度數(shù)值解的情況下，如何大幅度加速油藏模擬運(yùn)算速度是亟待解決的問(wèn)題.

模型降階技術(shù)最早出現(xiàn)在自動(dòng)控制和電路系統(tǒng)領(lǐng)域，它的任務(wù)是減少系統(tǒng)狀態(tài)空間向量的維數(shù)，同時(shí)保持系統(tǒng)的輸入、輸出特性.本征正交分解法（Proper Orthogonal Decomposition，簡(jiǎn)記為POD）是非線性系統(tǒng)模型降階中最為廣泛使用的方法.該方法所需的計(jì)算代價(jià)較低，具有一定穩(wěn)定性，同時(shí)又能保持原始系統(tǒng)的一些基本性質(zhì).目前POD方法已經(jīng)成功地應(yīng)用于計(jì)算流體力學(xué)[6]、結(jié)構(gòu)力學(xué)[7]、氣象學(xué)[8]及數(shù)字信號(hào)處理[9]等領(lǐng)域.本文將POD模型降階方法應(yīng)用于油藏模擬器，可大幅度降低油藏模型維數(shù)，從而減少計(jì)算時(shí)間，提高運(yùn)算速度.

1 油藏模型的控制方程

本文利用空間上的離散將油藏模型的數(shù)學(xué)方程轉(zhuǎn)換成狀態(tài)空間方程的形式，以便說(shuō)明POD方法的降階過(guò)程.采用二維油水兩相油藏模型，假設(shè)油、水兩相不發(fā)生物質(zhì)交換、過(guò)程等溫、流體微可壓縮，應(yīng)用質(zhì)量守恒方程和達(dá)西定律可得[10]

式中：K為滲透率張量；μ為流體粘度；kr為相對(duì)滲透率；p為壓力；g為重力加速度；d為深度；ρ為流體密度；φ為孔隙度；S為流體飽和度；t為時(shí)間；q000為源匯相，表示單位體積的流速；上標(biāo)i∈{o,w}分別表示油相和水相.等式（1）中有四個(gè)未知量po,pw,So,Sw，利用輔助方程（2）、（3），消去pw,So，使等式中只包含po,Sw兩個(gè)狀態(tài)變量.

式中：pc(Sw)為油水兩相的毛管壓力.

我們這里考慮較簡(jiǎn)單的情況，忽略重力和毛管力.利用五點(diǎn)塊中心有限差分格式在空間上進(jìn)行離散，可得非線性一階微分方程[11].

式中：向量p和s分別為網(wǎng)格中心油相壓力po和水相飽和度Sw；V為累積矩陣；T為傳導(dǎo)矩陣；F為分流量矩陣；向量qwell,t為井的油水總流量.

定義狀態(tài)向量x，輸入向量u，輸出向量y分別為

等式（4）可寫成狀態(tài)空間方程的形式[11]

在控制系統(tǒng)中，Ac稱為系統(tǒng)矩陣，Bc稱為輸入矩陣，C稱為輸出矩陣，D稱為直接傳遞矩陣.由于矩陣V，T，F(xiàn)中的元素都是狀態(tài)變量x的函數(shù)，則該系統(tǒng)為非線性系統(tǒng).

下面采用全隱式方法求解控制方程（8），用?x/?t近似導(dǎo)數(shù)dx/dt，則有

設(shè)xk表示tk時(shí)刻的狀態(tài)變量值，xk?1表示tk?1時(shí)刻的狀態(tài)變量值，則有

令fk=xk?1+?tfc(uk,xk)，可得

由于等式（12）左右兩端均含xk，則它為非線性方程組，一般采用牛頓迭代法求解.但當(dāng)變量xk的維數(shù)很高時(shí)，求解過(guò)程非常耗時(shí)，計(jì)算代價(jià)較大，為此，我們將POD模型降階方法應(yīng)用于油藏模擬，降低xk的維數(shù)，進(jìn)而提高運(yùn)算速度.

2 POD模型降階方法

POD方法也被稱為經(jīng)驗(yàn)正交函數(shù)法和Karhunen–Lo`eve分解法.該方法利用給定的理論或者實(shí)驗(yàn)得到的數(shù)據(jù)來(lái)構(gòu)造一組最優(yōu)正交基，達(dá)到降階目的.下面我們給出該方法的基本原理及在油藏?cái)?shù)值模擬中的具體應(yīng)用，關(guān)于POD方法更詳細(xì)的介紹可參考文獻(xiàn)[12-15].

2.1 POD方法的基本原理[16]

假設(shè)φ1,φ2,…,φn是Rn中數(shù)據(jù)集{x1,x2,…,xm}的標(biāo)準(zhǔn)正交基向量，則對(duì)任意樣本向量xi∈Rn，就可以由基底按如下形式線性表出

其中cji=.若選擇前k個(gè)基向量來(lái)逼近樣本向量xi，則有

對(duì)于獲得的數(shù)據(jù)集{x1,x2,…,xm}，POD方法的目標(biāo)是找到一組最優(yōu)基向量{φ1,φ2,…,φn}，以使誤差ε2達(dá)到最小.引入Lagrange因子uij(i,j=k+1,k+2,…,n)，并且構(gòu)造Lagrange函數(shù)

（16）式兩端對(duì)φj求偏導(dǎo)，得到

其中Un?k=[uk+1uk+2…un].為求最優(yōu)解，讓（18）式左端為0，則有

其中Λ是對(duì)角矩陣.對(duì)（19）式兩端右乘矩陣P，并結(jié)合（20）式，就有

由（21）式可知，矩陣Λ的對(duì)角元素由矩陣XXT的特征值λi構(gòu)成，矩陣Φn?kP由矩陣XXT相應(yīng)于λi的特征向量構(gòu)成.實(shí)際上，矩陣Λ的對(duì)角元素是矩陣XT的奇異值σi的平方.考慮到正交矩陣在Frobenius范數(shù)k?kF下是保范的，因此得到

為了獲得最小誤差，矩陣Λ的對(duì)角元素需要選擇矩陣XT的最后n?k個(gè)奇異值.這樣，就有minε2=

由以上討論可知，當(dāng)選擇矩陣XT的右奇異值向優(yōu)性條件，而且也獲得了最小誤差，其最小誤差為矩陣XT的最后n?k個(gè)奇異值的平方和.

通常在輸入、輸出系統(tǒng)中，較大的特征值對(duì)應(yīng)系統(tǒng)的主要特征.為此，我們選擇的k個(gè)POD基向要能代表原始向量較多的特征相關(guān)矩陣按降序排列的特征值.

假如I(k)≥d%，則k個(gè)POD基向量保持原始數(shù)據(jù)集的d%特征信息，其中k滿足

2.2 POD方法在油藏模擬中的應(yīng)用

為了構(gòu)造POD基向量，首先要運(yùn)行全階油藏模擬器（也稱為訓(xùn)練過(guò)程），并保存每個(gè)時(shí)間步的狀態(tài)向量x（也稱為快照，包括所有網(wǎng)格的油相壓力po和水相飽和度Sw），由于壓力和飽和度具有不同的物理性質(zhì)，我們分別用矩陣Xp，XS保存po，Sw（以下簡(jiǎn)記為p和S）.

假設(shè)油藏模型的網(wǎng)格數(shù)為N，則矩陣Xp和XS中的每個(gè)向量（上標(biāo)i表示快照的數(shù)目）都是N維，而系統(tǒng)的狀態(tài)向量x的維數(shù)為n=2N.快照獲得后，需要計(jì)算快照的均值

并且數(shù)據(jù)矩陣Xp和XS中的每個(gè)快照都減去均值

將式（27）去掉下標(biāo)l代入（12）中，可得

等式兩端同時(shí)左乘ΦT

這就實(shí)現(xiàn)了對(duì)原系統(tǒng)降階的目的，狀態(tài)變量的維數(shù)從n=2N減少為l=lp+ls.再采用牛頓迭代法求解非線性方程組（29），由于狀態(tài)變量的維數(shù)較低，可大大減少計(jì)算量.

3 算例驗(yàn)證

本文采用文獻(xiàn)[11]中的算例.該算例描述的是一個(gè)二維油水兩相各向異性的方形油藏，網(wǎng)格劃分為21×21，滲透率和孔隙度的分布如圖1、2所示.油藏模型相關(guān)參數(shù)：油層厚度h=2 m，網(wǎng)格的長(zhǎng)、寬?x=?y=33.33m，原油粘度μo=5.0×10?4Pa·s，地層水粘度μw=1.0×10?3Pa·s，綜合壓縮系數(shù)ct=3.0×10?9Pa?1，原始地層壓力pi=30×106Pa，井孔半徑rwell=0.114 m，油相端點(diǎn)相對(duì)滲透率=0.9，水相端點(diǎn)相對(duì)滲透率=0.6，油相Corey指數(shù)no=2.0，水相Corey指數(shù)nw=2.0，殘余油飽和度Sor=0.2，束縛水飽和度Swc=0.2.采用反五點(diǎn)法井網(wǎng)生產(chǎn)，中心一口注水井，四角為四口生產(chǎn)井，忽略重力和毛管力影響.

圖1 油藏模型的滲透率分布

圖2 油藏模型的孔隙度分布

我們采用Matlab編輯的全隱式處理的數(shù)模程序模擬該算例，修改源代碼實(shí)現(xiàn)POD模型降階過(guò)程，驗(yàn)證該方法的有效性.

3.1 訓(xùn)練過(guò)程

生產(chǎn)井給定井底壓力為295×105Pa，注入井給定井底流量為0.0015m3/s，運(yùn)行全階模擬器模擬1400天，保存66個(gè)時(shí)間步的結(jié)果.壓力矩陣保留37個(gè)奇異值，飽和度矩陣保留38個(gè)奇異值，最后獲得POD的基矩陣Φl的維數(shù)為2N×l，其中l(wèi)=37+38=75.這就意味著全階模擬器要求解2N=882個(gè)未知變量，而降階后只需求解75個(gè)變量.

訓(xùn)練過(guò)程中，全階油藏模擬器與應(yīng)用POD方法的降階模擬器計(jì)算結(jié)果對(duì)比見(jiàn)圖3，4.

在文中我們采用平均相對(duì)誤差來(lái)衡量近似值的精度.如每口生產(chǎn)井產(chǎn)油量的平均相對(duì)誤差定義為

式中：i表示時(shí)間步；nt表示時(shí)間步總數(shù)；表示生產(chǎn)井m第i時(shí)間步全階模擬器的產(chǎn)油量；表示生產(chǎn)井m第i時(shí)間步POD降階模擬器的產(chǎn)油量.同理也可定義每口生產(chǎn)井產(chǎn)水量的平均相對(duì)誤差.

圖3 四口生產(chǎn)井的產(chǎn)油量對(duì)比（訓(xùn)練過(guò)程）

圖4 四口生產(chǎn)井的產(chǎn)水量對(duì)比（訓(xùn)練過(guò)程）

訓(xùn)練過(guò)程，四口生產(chǎn)井產(chǎn)油量、產(chǎn)水量的平均相對(duì)誤差如表1.

表1 四口生產(chǎn)井產(chǎn)油量、產(chǎn)水量的平均相對(duì)誤差（訓(xùn)練過(guò)程）

以上結(jié)果顯示，在訓(xùn)練過(guò)程中，降階模擬器與全階模擬器四口生產(chǎn)井的產(chǎn)油量和產(chǎn)水量基本完全相同，產(chǎn)生的平均相對(duì)誤差都非常小，但模擬時(shí)間提高了近3倍，全階模擬器的運(yùn)行時(shí)間為35.527 s，而降階模擬器的運(yùn)行時(shí)間為12.019 s.

3.2 預(yù)測(cè)過(guò)程

下面我們用訓(xùn)練過(guò)程獲得的基矩陣Φl來(lái)驗(yàn)證POD降階模型的預(yù)測(cè)能力.此時(shí)生產(chǎn)井的井底壓力改為285×105Pa，注入井的井底流量保持不變，全階模擬器與降階模擬器的計(jì)算結(jié)果對(duì)比見(jiàn)圖5、6.

圖5 四口生產(chǎn)井的產(chǎn)油量對(duì)比（預(yù)測(cè)過(guò)程）

圖6 四口生產(chǎn)井的產(chǎn)水量對(duì)比（預(yù)測(cè)過(guò)程）

預(yù)測(cè)過(guò)程，四口生產(chǎn)井產(chǎn)油量、產(chǎn)水量的平均相對(duì)誤差見(jiàn)表2.

表2 四口生產(chǎn)井產(chǎn)油量、產(chǎn)水量的平均相對(duì)誤差（預(yù)測(cè)過(guò)程）

上面結(jié)果顯示，當(dāng)預(yù)測(cè)過(guò)程生產(chǎn)制度與訓(xùn)練過(guò)程生產(chǎn)制度不同時(shí)，降階與全階模擬器輸出結(jié)果的平均相對(duì)誤差有所提高，但仍控制在5%的合理范圍內(nèi).此時(shí)，模擬時(shí)間也提高了將近3倍，全階模擬器的運(yùn)行時(shí)間為35.851 s，降階模擬器的運(yùn)行時(shí)間為12.404 s.

4 結(jié)論

（1）將油水兩相滲流的偏微分方程，利用五點(diǎn)塊中心有限差分格式在空間上進(jìn)行離散，轉(zhuǎn)化成狀態(tài)空間方程的形式；

（2）POD模型降階方法應(yīng)用于油藏模擬器，可大幅度減少油藏模型維數(shù)，實(shí)例顯示在保證一定精度的條件下，模擬器的運(yùn)算速度提高近3倍，驗(yàn)證了該方法的有效性；

（3）當(dāng)訓(xùn)練與檢測(cè)過(guò)程的生產(chǎn)制度不同時(shí)，降階模擬器的計(jì)算結(jié)果平均相對(duì)誤差有所提高，并且兩個(gè)過(guò)程的生產(chǎn)制度差別越大，產(chǎn)生誤差也越大，所以POD模型降階方法還需進(jìn)一步改進(jìn).