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

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

0 引言

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

目前，中國新開發(fā)油田中低滲、特低滲儲量比例高（約占60%～70%），已開發(fā)主力油田進(jìn)入高含水、高采出程度的開發(fā)后期[6-7]，油藏開發(fā)難度加大。為了進(jìn)一步獲得經(jīng)濟(jì)有效的產(chǎn)量，需深化油藏描述，建立更精細(xì)的油藏?cái)?shù)值模型[8-9]，這些模型比以往的地質(zhì)模型網(wǎng)格數(shù)更多、生產(chǎn)歷史更長、井的數(shù)目和措施更多、油氣分布更為復(fù)雜，極大地增加了油藏?cái)?shù)值模擬歷史擬合和調(diào)整方案開發(fā)趨勢預(yù)測的工作量、難度并更耗時(shí)[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——時(shí)間，s；pw，po，pg——水相、油相、氣相壓力，MPa；g——重力加速度，m/s2；Z——深度，m；pcow，pcgo——油水、氣油界面毛管壓力，MPa。

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

其中

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

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

2 塊壓縮存儲格式

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

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

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

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

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

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

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

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

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

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

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

有效節(jié)點(diǎn)壓縮存儲雖然能有效減少內(nèi)存的損耗，但壓縮后的稀疏矩陣中仍含有大量的零元素?？紤]到多變量黑油模型的特點(diǎn)，本文以塊為單位存儲系數(shù)矩陣中的非零元素，根據(jù)矩陣中非零塊元素（塊中有一個(gè)非零元素即認(rèn)為是非零塊元素）的位置對稱性，只記錄下三角矩陣非零塊元素的地址。本文采用 3個(gè)雙精度指針數(shù)組 E、T及 T?分別存儲系數(shù)矩陣中的對角線塊元素、下三角塊元素及上三角塊元素，并利用 3個(gè)整型數(shù)組記錄非零塊元素的地址，實(shí)現(xiàn)稀疏矩陣塊壓縮排列：①Nco（SI），下三角矩陣中每行非零塊元素的總數(shù)（SI為塊行號）；②Nod（SI）：下三角矩陣中每塊行第1個(gè)非零塊元素的位置序號；③Icn（SR）：下三角矩陣中非零塊元素對應(yīng)的列序號（SR為非零塊元素的累計(jì)計(jì)數(shù)序號）。經(jīng)有效節(jié)點(diǎn)壓縮后，黑油模型（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個(gè)數(shù)組：1個(gè)實(shí)數(shù)型數(shù)組存儲稀疏矩陣的非零元素，2個(gè)整型數(shù)組分別記錄相應(yīng)非零元素行和列的位置，若 n階稀疏矩陣中非零塊元素個(gè)數(shù)為 nt，其中每個(gè)非零塊元素均為3階子矩陣，則該壓縮格式所需存儲空間為：9ntLd+(3n+1+9nt)Lint，其中 Ld為一個(gè)雙精度數(shù)長度（8 B），Lint為一個(gè)整數(shù)的長度（2 B）；對于相同的n階nt個(gè)非零 3階塊元素的稀疏矩陣，按本文提出的塊壓縮存儲需要存儲空間：9ntLd+ [(3n + nt)/2+1] Lint；若不利用壓縮存儲格式需要的存儲空間是n2Ld。以三相黑油模型在不同有效網(wǎng)格上形成的系數(shù)矩陣為例，塊壓縮存儲比行壓縮存儲耗費(fèi)內(nèi)存少，且矩陣規(guī)模越大節(jié)約的內(nèi)存越多（見表5）。

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

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

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

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

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

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

①令m=1；

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

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

③對3階子陣T及T ?中的元素進(jìn)行與第1步類似的計(jì)算，得到相對應(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ù)處理迭代法的一般流程如下。

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

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

③計(jì)算 wj=M?1Avj；

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

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

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

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

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

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

國際黑油模型標(biāo)準(zhǔn)考題SPE10描述了大規(guī)模嚴(yán)重非均質(zhì)水驅(qū)油藏的生產(chǎn)動(dòng)態(tài)，是反映油藏?cái)?shù)值模擬技術(shù)的規(guī)?；⒂?jì)算精度和計(jì)算速度的重要考題。該算例為油水兩相模型，油藏埋深為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ā)動(dòng)態(tài)，總節(jié)點(diǎn)數(shù)為1 122 000。在處理器為Intel Xeon X5660 2.80GHz的計(jì)算機(jī)上用塊GMRES預(yù)處理算法、商業(yè)模擬器A及B模擬SPE10算例的日產(chǎn)油、綜合含水和地層壓力等生產(chǎn)指標(biāo)曲線（見圖 3），可見，該算法的計(jì)算精度與商業(yè)油藏?cái)?shù)值模擬器基本一致。在該計(jì)算機(jī)上，塊GMRES預(yù)處理算法、商業(yè)模擬器A及B模擬SPE10算例的總耗時(shí)分別為21.5 h、24.3 h和101.3 h，計(jì)算速度較快。說明塊 GMRES預(yù)處理算法能夠快速實(shí)現(xiàn)大規(guī)模的油藏?cái)?shù)值模擬，結(jié)果可靠，具有較強(qiáng)的實(shí)用性。

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

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

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

5 結(jié)論

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

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

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

[1] 韓大匡, 陳欽雷, 閆存章. 油藏?cái)?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. 油藏?cái)?shù)值模擬自動(dòng)歷史擬合方法: 以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ì)算[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í)層析成像中的大型稀疏系數(shù)矩陣壓縮存儲和求解[J]. 地球物理學(xué)進(jìn)展, 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] 林小兵. 二維多相油藏?cái)?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.