王寶華，吳淑紅，韓大匡，桓冠仁，李巧云，李小波，李華，周久寧

（1. 中國石油勘探開發(fā)研究院；2. 提高石油采收率國家重點實驗室；3. 中國石油大學(xué)（北京））

0 引言

油藏數(shù)值模擬利用已知油藏地質(zhì)參數(shù)及動態(tài)開發(fā)數(shù)據(jù)再現(xiàn)油藏開發(fā)歷史、預(yù)測油藏開發(fā)趨勢、制定或調(diào)整開發(fā)方案[1-3]，已成為油田開發(fā)研究必不可少的技術(shù)手段。油藏數(shù)值模擬的數(shù)學(xué)模型由多個非線性偏微分方程耦合而成，具有強非線性、強間斷性、強耦合性，求解這些微分方程組的時間通常占整個模擬時間的 70%～80%，隨著模型規(guī)模和復(fù)雜程度的增加，該比例還會增加，如何高效求解這些微分方程組成為數(shù)值模擬的核心問題[4-5]。

目前，中國新開發(fā)油田中低滲、特低滲儲量比例高（約占60%～70%），已開發(fā)主力油田進入高含水、高采出程度的開發(fā)后期[6-7]，油藏開發(fā)難度加大。為了進一步獲得經(jīng)濟有效的產(chǎn)量，需深化油藏描述，建立更精細(xì)的油藏數(shù)值模型[8-9]，這些模型比以往的地質(zhì)模型網(wǎng)格數(shù)更多、生產(chǎn)歷史更長、井的數(shù)目和措施更多、油氣分布更為復(fù)雜，極大地增加了油藏數(shù)值模擬歷史擬合和調(diào)整方案開發(fā)趨勢預(yù)測的工作量、難度并更耗時[10]。為了提高數(shù)值模擬效率，本文針對油藏模擬中最常用的黑油模型，研究大規(guī)模稀疏系數(shù)矩陣的壓縮存儲及求解方法。

1 基礎(chǔ)模型

黑油模型中油、氣、水3相基本微分方程為：

各相壓力和毛細(xì)管壓力有關(guān)，定義為：

式中 K——滲透率，10?3μm2；φ——孔隙度，%；Krw，Kro，Krg——水相、油相、氣相相對滲透率，f；μw，μo，μg——水相、油相、氣相黏度，mPa·s；ρw，ρo，ρg——水相、油相、氣相密度，g/cm3；ρo,o——油相中油組分的密度，g/cm3；ρg,o——油相中氣組分的密度，g/cm3；Sw，So，Sg——水相、油相、氣相飽和度，%；t——時間，s；pw，po，pg——水相、油相、氣相壓力，MPa；g——重力加速度，m/s2；Z——深度，m；pcow，pcgo——油水、氣油界面毛管壓力，MPa。

假設(shè)三維模型橫向上有nx個網(wǎng)格，縱向上ny個網(wǎng)格（每個平面上有nxy個網(wǎng)格，nxy=nxny），垂向上有nz個網(wǎng)格，模型總網(wǎng)格數(shù)為n（n=nxnynz），在按行自然排序的塊中心網(wǎng)格上，利用控制體積全隱式求解格式離散黑油模型（1）式—（3）式，得到 n階帶狀稀疏系數(shù)矩陣：

其中

式中 B，C，D——組分差分方程中水飽和度、油飽和度和油相壓力 3個變量在相關(guān)有效節(jié)點上的系數(shù)；下標(biāo)：w，o，g——水、油、氣組分。

根據(jù)黑油模型的數(shù)學(xué)特性，將具有拋物方程性質(zhì)的油相壓力方程和具有雙曲型方程性質(zhì)的飽和度方程結(jié)合在一起，以塊形式進行運算能夠避免小擾動對模擬問題的干擾，增強穩(wěn)定性，減少迭代次數(shù)，提高模擬速度。隨著油藏模擬規(guī)模的增加，稀疏系數(shù)矩陣中零元素大量增加，遠(yuǎn)超過非零元素的數(shù)目，稀疏矩陣存儲方式?jīng)Q定了內(nèi)存的消耗量，直接影響運算時間。

2 塊壓縮存儲格式

油藏數(shù)值模擬中采用的壓縮存儲格式有很多，如點的坐標(biāo)壓縮存儲、行壓縮存儲、修正行壓縮存儲、對角線壓縮存儲及由點的行壓縮存儲演化出的塊行壓縮存儲等[2,11]。本文主要研究塊的壓縮存儲格式，主要由有效節(jié)點壓縮存儲和塊壓縮存儲兩部分構(gòu)成，與行壓縮存儲格式相比，塊壓縮存儲格式能夠節(jié)約更多的內(nèi)存，減少非零元素的搜索次數(shù)。

2.1 有效節(jié)點壓縮存儲技術(shù)

在實際油藏模型中，地層中存在大量死結(jié)點（如孔隙度為零），據(jù)此特點，在稀疏矩陣壓縮存儲時，只存儲有效節(jié)點上的相應(yīng)信息從而達到矩陣降階、節(jié)約內(nèi)存的目的。以一個3行3列的二維網(wǎng)格為例，其中第5、6、8、9這4個節(jié)點為有效節(jié)點，其余為死節(jié)點（見圖1a），若全部存儲，三變量黑油模型的系數(shù)矩陣將達到27階；若將這4個有效節(jié)點重新排序（見圖1b），則黑油模型系數(shù)矩陣降至12階。

為將此系數(shù)矩陣降階，定義壓縮排列數(shù)組：

圖1 自然排序（a）和有效節(jié)點排序（b）的二維網(wǎng)格

其中 S =（J?1）Nx+I

式中 N——壓縮排列序號，1≤N≤ne；ne——有效節(jié)點總數(shù)；Nc——有效網(wǎng)格節(jié)點映射；S——常規(guī)排列序號；J——縱向網(wǎng)格序號，1≤J≤ny；I——橫向網(wǎng)格序號，1≤I≤nx；nx，ny——橫、縱方向的網(wǎng)格總數(shù)。

利用（6）式將常規(guī)排列中的死節(jié)點去除并將有效節(jié)點重新排序（見表 1），死節(jié)點置 0，有效節(jié)點按自然順序排列。

表1 有效節(jié)點排序表

利用上述方法壓縮有效節(jié)點能夠大幅降低稀疏系數(shù)矩陣的階數(shù)，節(jié)約模擬內(nèi)存，如中國東部某陸相非均質(zhì)油田數(shù)值模型中，總結(jié)點數(shù)約 40×104，有效節(jié)點數(shù)約14×104，若全部存儲需要0.200 GB的內(nèi)存，若用本文壓縮排列方法只需 0.065 GB的內(nèi)存，節(jié)約了近70%的內(nèi)存空間。

2.2 塊壓縮存儲技術(shù)

有效節(jié)點壓縮存儲雖然能有效減少內(nèi)存的損耗，但壓縮后的稀疏矩陣中仍含有大量的零元素?？紤]到多變量黑油模型的特點，本文以塊為單位存儲系數(shù)矩陣中的非零元素，根據(jù)矩陣中非零塊元素（塊中有一個非零元素即認(rèn)為是非零塊元素）的位置對稱性，只記錄下三角矩陣非零塊元素的地址。本文采用 3個雙精度指針數(shù)組 E、T及 T?分別存儲系數(shù)矩陣中的對角線塊元素、下三角塊元素及上三角塊元素，并利用 3個整型數(shù)組記錄非零塊元素的地址，實現(xiàn)稀疏矩陣塊壓縮排列：①Nco（SI），下三角矩陣中每行非零塊元素的總數(shù)（SI為塊行號）；②Nod（SI）：下三角矩陣中每塊行第1個非零塊元素的位置序號；③Icn（SR）：下三角矩陣中非零塊元素對應(yīng)的列序號（SR為非零塊元素的累計計數(shù)序號）。經(jīng)有效節(jié)點壓縮后，黑油模型（1）式—（3）式離散形成如下形式的矩陣（以 4階矩陣為例）：

A中各元素為如（5）式所示的3階子陣。利用塊壓縮存儲格式，矩陣 A的存儲見表 2，各元素按行依次存儲。

表2 矩陣A的塊壓縮存儲表

非零塊元素的地址由數(shù)組Nco、Nod及Icn記錄（對應(yīng)于下三角壓縮存儲數(shù)組T），見表3、表4。

表3 塊壓縮存儲的Nco和Nod數(shù)組

表4 塊壓縮存儲的Icn數(shù)組

2.3 不同壓縮存儲方法存儲量及搜索次數(shù)對比

2.3.1 存儲量對比

不同壓縮存儲格式存儲稀疏矩陣非零元素所需的數(shù)組長度大致相同，但記錄非零元素地址的數(shù)組差別很大，如行壓縮存儲格式存儲稀疏矩陣需要3個數(shù)組：1個實數(shù)型數(shù)組存儲稀疏矩陣的非零元素，2個整型數(shù)組分別記錄相應(yīng)非零元素行和列的位置，若 n階稀疏矩陣中非零塊元素個數(shù)為 nt，其中每個非零塊元素均為3階子矩陣，則該壓縮格式所需存儲空間為：9ntLd+(3n+1+9nt)Lint，其中 Ld為一個雙精度數(shù)長度（8 B），Lint為一個整數(shù)的長度（2 B）；對于相同的n階nt個非零 3階塊元素的稀疏矩陣，按本文提出的塊壓縮存儲需要存儲空間：9ntLd+ [(3n + nt)/2+1] Lint；若不利用壓縮存儲格式需要的存儲空間是n2Ld。以三相黑油模型在不同有效網(wǎng)格上形成的系數(shù)矩陣為例，塊壓縮存儲比行壓縮存儲耗費內(nèi)存少，且矩陣規(guī)模越大節(jié)約的內(nèi)存越多（見表5）。

表5 塊壓縮存儲、行壓縮存儲及非壓縮存儲占用的內(nèi)存

2.3.2 搜索次數(shù)的對比

在模擬中，對于非零元素的搜索，塊壓縮存儲遠(yuǎn)遠(yuǎn)少于行壓縮存儲，以n階nt個非零3階塊元素的稀疏矩陣為例，在行壓縮格式中需要搜索 9nt個非零元素，而在塊壓縮格式中只需搜索下三角矩陣中非零塊元素，即只需搜索(nt?n)/2個非零塊元素，行壓縮存儲需要的搜索次數(shù)大約是塊壓縮存儲的 18nt/(nt－n)倍。仍以三相黑油模型在不同有效網(wǎng)格上離散的系數(shù)矩陣為例，塊壓縮存儲需要搜索的非零元素個數(shù)遠(yuǎn)少于行壓縮存儲，大約只有行壓縮存儲的5%（見表6）。

表6 塊壓縮存儲、行壓縮存儲非零元素搜索次數(shù)

3 塊廣義極小余量 GMRES預(yù)處理迭代法

三維三相黑油模型（1）式—（3）式離散后形成大規(guī)模稀疏線性方程組，求解這些線性方程組大約占整體模擬時間的80%左右，而且隨著模擬規(guī)模的增加，該比例還會增大。目前求解大型稀疏矩陣較為成熟的方法是預(yù)處理迭代方法，如不完全LU分解（incomplete LU factorization）預(yù)處理廣義極小殘量迭代法（generalized minimal residual method，GMRES），該算法利用 Krylov子空間矢量的最小殘量進行迭代求解，具有收斂速度快、穩(wěn)定性好等優(yōu)點[12-13]。由于油藏問題中的系數(shù)矩陣一般很難滿足點的對角占優(yōu)性，本文考慮以節(jié)點為單元進行ILU分解GMRES迭代。

3.1 塊不完全LU分解

塊不完全LU分解就是預(yù)定分解矩陣L、U的非零指標(biāo)集G，在塊LU分解中置G以外子塊為零，形成系數(shù)矩陣A的近似分解：A（G）=L（G）U（G）[14]。前文中已將形如（4）式的帶狀矩陣 A = (aij)n×n的主對角、下三角及上三角非零塊元素分別存入E、T、及T ?數(shù)組中，利用這3個數(shù)組對A進行塊LU分解時只有A的非零元素參與計算而且不產(chǎn)生新的非零元素，如此計算得到的L、U矩陣為A的塊不完全LU分解，具體計算步驟如下。

①令m=1；

②將對角線上第m個塊矩陣進行標(biāo)準(zhǔn)化處理，分解為上三角塊與下三角塊的乘積，下三角塊中的元素lji與上三角塊中的元素uij由下式計算：

其中，Cij為塊對角矩陣的元素。

③對3階子陣T及T ?中的元素進行與第1步類似的計算，得到相對應(yīng)的L、U矩陣中的子陣的值；

④令m=m+1，返回②直至m=n。

3.2 塊廣義極小殘量迭代（GMRES）法

設(shè) n階線性方程組 Ax=b（A如（4）式，其元素 aij如（5）式，x為未知量，其初始近似解向量為x0，第l次迭代的近似解向量為xl，相應(yīng)的殘差為rl=b?Axl，l=0，1，…。GMRES算法是通過求使 rl極小的xl∈κl來逼近Ax=b的精確解，其中κl為Krylov子空間：為首次殘差，設(shè)Vl=(v1，v2，…，vl)為κl的標(biāo)準(zhǔn)正交基，且AVl=Vl+1Hl，其中Hl為上Hessenberg矩陣，xl可寫成xl=Vlyl（yl為l維向量），則rl的極小值問題轉(zhuǎn)化為極小問題：

其中02||||β=r，e1=(1，0，0，…，0)?[13]。這種基于Krylov子空間的迭代法收斂速度快、精度高且穩(wěn)定性好，塊GMRES預(yù)處理迭代法的一般流程如下。

①計算r0=M?1（b?Ax0），此處M=LU是A的ILU分解陣，以及標(biāo)準(zhǔn)正交基第1個元素10β=vr，其中矩陣每個塊aij與向量塊xj的乘積如下：

②循環(huán) j=1，2，…，l；

③計算 wj=M?1Avj；

④對于i=1，2，…，j，計算Hessenberg矩陣中的向量 hij=（wj, vi）；

⑤計算Hessenberg矩陣中的向量hj+1,j= | |wj||2和vj+1=wj||wj||2；

⑥更新xl=x0+Vlyl，其中yl滿足（7）式；

⑦如果收斂，停止循環(huán)；否則，令x0= xl，返回①。

步驟②—⑥中矩陣與向量的乘積及向量與向量的運算均為塊矩陣之間的運算。

4 大規(guī)模水驅(qū)油藏算例模擬

國際黑油模型標(biāo)準(zhǔn)考題SPE10描述了大規(guī)模嚴(yán)重非均質(zhì)水驅(qū)油藏的生產(chǎn)動態(tài)，是反映油藏數(shù)值模擬技術(shù)的規(guī)?；?、計算精度和計算速度的重要考題。該算例為油水兩相模型，油藏埋深為3 657.6 m，初始油藏壓力41.37 MPa，地面油密度、氣體密度及地面水密度分別為 0.849 g/cm3、1.000×10?4g/cm3及 1.024 g/cm3，地層水及地下原油體積系數(shù)均為 1.01，水和巖石壓縮系數(shù)分別為 4.41×10?5MPa?1和 1.47×10?5MPa?1，水及地下原油黏度分別為0.3 mPa·s和3.0 mPa·s。

在 SPE10算例中滲透率的最大值為 20 000×10?3μm2，最小值為 0.000 66×10?3μm2，平均值為364.52×10?3μm2，油層滲透率極差為 23.3，垂向與平面滲透率的比值變化較大，其中，河道為0.3，基底的為 0.001。該算例中孔隙度最大值為 0.5，平均值為0.174 9，部分單層孔隙度見圖2。

圖2 SPE10算例部分單層孔隙度分布

SPE10模型模擬的是366 m×671 m×51.85 m范圍內(nèi)1口注水井和4口生產(chǎn)井的開發(fā)動態(tài)，總節(jié)點數(shù)為1 122 000。在處理器為Intel Xeon X5660 2.80GHz的計算機上用塊GMRES預(yù)處理算法、商業(yè)模擬器A及B模擬SPE10算例的日產(chǎn)油、綜合含水和地層壓力等生產(chǎn)指標(biāo)曲線（見圖 3），可見，該算法的計算精度與商業(yè)油藏數(shù)值模擬器基本一致。在該計算機上，塊GMRES預(yù)處理算法、商業(yè)模擬器A及B模擬SPE10算例的總耗時分別為21.5 h、24.3 h和101.3 h，計算速度較快。說明塊 GMRES預(yù)處理算法能夠快速實現(xiàn)大規(guī)模的油藏數(shù)值模擬，結(jié)果可靠，具有較強的實用性。

圖3 SPE10全油田日產(chǎn)油、綜合含水及地層壓力曲線

圖4 為采用塊GMRES預(yù)處理算法模擬SPE10算例得到的不同時間段注采井間水線推進剖面。該方法可以精確模擬注水過程中水線推進到油藏并逐步到達油井的過程，對于嚴(yán)重非均質(zhì)油藏，可以清晰地分析描述注入水在不同滲透層推進的速度以及沿高滲透層和高滲區(qū)域的指進過程。

圖4 SPE10算例在不同時刻含水飽和度分布（P1、P3為生產(chǎn)井，I1為注水井）

5 結(jié)論

針對三相黑油模型模擬時形成的大規(guī)模稀疏系數(shù)矩陣，提出了有效節(jié)點壓縮和塊壓縮存儲相結(jié)合的存儲格式，分別采用3個實數(shù)組和3個整型數(shù)組存儲稀疏矩陣的非零塊元素及其地址信息，減少了內(nèi)存的占用并降低了計算時的搜索次數(shù)。

在求解黑油模型的稀疏線性方程組時，先用塊ILU分解改善系數(shù)矩陣的特征值分布，得到性質(zhì)更好的矩陣，再用塊 GMRES算法求解改良的線性方程組。由于黑油模型的系數(shù)矩陣具有塊對角占優(yōu)性質(zhì)，塊GMRES預(yù)處理法具有良好的穩(wěn)定性和收斂性，能夠快速、精確實現(xiàn)大規(guī)模油藏數(shù)值模擬。

塊壓縮存儲與塊 ILU GMRES算法相結(jié)合可模擬大規(guī)模油藏的生產(chǎn)動態(tài)，計算結(jié)果與商業(yè)軟件基本一致，速度較快且耗費內(nèi)存較少，是一種高效、可靠的求解技術(shù)。

[1] 韓大匡, 陳欽雷, 閆存章. 油藏數(shù)值模擬基礎(chǔ)[M]. 北京: 石油工業(yè)出版社, 1999: 34-44, 241-242.Han Dakuang, Chen Qinlei, Yan Cunzhang. The basis of reservoir simulation[M]. Beijing: Petroleum Industry Press, 1999: 34-44, 241-242.

[2] Chen Zhangxin, Huan Guanren, Ma Yuanle. Computational methods for multiphase flows in porous media[M]. Philadelphia: Society for Industrial and Applied Mathematics, 2006: 1-2.

[3] Kazemi A, Stephen K D. 油藏數(shù)值模擬自動歷史擬合方法: 以Nelson油田為例[J]. 石油勘探與開發(fā), 2012, 39(3): 326-337.Kazemi A, Stephen K D. Schemes for automatic history matching of reservoir modeling: A case of Nelson oilfield in the UK[J].Petroleum Exploration and Development, 2012, 39(3): 326-337.

[4] Li Xiaobo, Wu Shuhong, Song Jie, et al. Numerical simulation of pore- scale flow in chemical flooding process[J]. Theor. Appl. Mech.Lett., 2011, 1(2): 022008.

[5] Wang Baohua, Lü Shujuan, Lu Qishao, et al. Long time behavior of finite difference solution for a semilinear parabolic equation[J].DCDIS-A, 2007, 14(S2): 141-147.

[6] 紀(jì)淑紅, 田昌炳, 石成方, 等. 高含水階段重新認(rèn)識水驅(qū)油效率[J]. 石油勘探與開發(fā), 2012, 39(3): 338-345.Ji Shuhong, Tian Changbing, Shi Chengfang, et al. New understanding on water-oil displacement efficiency in a high water-cut stage[J]. Petroleum Exploration and Development, 2012,39(3): 338-345.

[7] 李巧云, 張吉群, 鄧寶榮, 等. 高含水油田層系重組方案的灰色決策優(yōu)選法[J]. 石油勘探與開發(fā), 2011, 38(4): 463-468.Li Qiaoyun, Zhang Jiqun, Deng Baorong, et al. Grey decision-making theory in the optimization of strata series recombination programs of high water-cut oilfields[J]. Petroleum Exploration and Development, 2011, 38(4): 463-468.

[8] Wu Shuhong, Han Min, Ma Desheng, et al. Streamflooding to enhance recovery of a waterflooded light-oil reservoir[J]. JPT,2012(1): 64-66.

[9] 韓大匡. 關(guān)于高含水油田二次開發(fā)理念、對策和技術(shù)路線的探討[J]. 石油勘探與開發(fā), 2010, 37(5)：583-591.Han Dakuang. Discussions on concepts, countermeasures and technical routes for the redevelopment of high water-cut oilfields[J].Petroleum Exploration and Development, 2010, 37(5): 583-591.

[10] 鄒存友, 常毓文, 王國輝, 等. 水驅(qū)開發(fā)油田合理油水井?dāng)?shù)比的計算[J]. 石油勘探與開發(fā), 2011, 38(2): 211-215.Zou Cunyou, Chang Yuwen, Wang Guohui, et al. Calculation on a reasonable production-injection well ratio in waterflooding oilfields[J].Petroleum Exploration and Development, 2011, 38(2): 211-215.

[11] 成谷, 張寶金. 反射地震走時層析成像中的大型稀疏系數(shù)矩陣壓縮存儲和求解[J]. 地球物理學(xué)進展, 2008, 23(3): 674-680.Cheng Gu, Zhang Baojin. Compression storage and solution of large sparse matrix in traveltime tomography of reflection seismic data[J].Progress in Geophysics, 2008, 23(3): 674-680.

[12] Saad Y, Gmres M S. A generalized minimal residual algorithm for solving nonsymmetric linear systems[J]. Siam J. Sci. Stat. Comput.,1986, 7(3): 856-869.

[13] Li Wenjun, Chen Zhangxi, Richard E E, et al. Comparision of the GMRES and ORTHOMIN for the black oil model in porous media[J].Int. J. Numer. Meth. Fluids., 2005, 48(5): 501-519.

[14] 林小兵. 二維多相油藏數(shù)值模擬的 BILUCG算法[J]. 石油勘探與開發(fā), 1987, 14(1): 76-84.Lin Xiaobing. BILUCG algorithm for 2-D multiphase reservoir simulation[J]. Petroleum Exploration and Development, 1987, 14(1):76-84.