朱舟元，李明輝，雷征東，陳掌星

1 中國石油大學(xué)(北京)石油工程學(xué)院，北京 102200 2 中國石油大學(xué)(北京)非常規(guī)天然氣研究院，北京 102200 3 中國石油勘探開發(fā)研究院，北京 100083 4 中國石油大學(xué)(北京)油氣資源與工程國家重點(diǎn)實(shí)驗(yàn)室，北京 102200

0 引言

復(fù)雜油氣藏采收過程的數(shù)值模擬往往包含了復(fù)雜的物理化學(xué)過程，油藏網(wǎng)格作為模擬多相多組分流動(dòng)的基本單元，它的尺寸大小對(duì)于復(fù)雜開采過程的模擬精度有著重要的影響。一般認(rèn)為網(wǎng)格尺寸越小，模擬結(jié)果的精度越高，越趨近于收斂的解。網(wǎng)格尺寸越大，數(shù)值誤差越大，造成數(shù)值彌散、耗散等現(xiàn)象，模擬結(jié)果越不精確[1]。故在可行的情況下，模擬必須采用足夠精細(xì)的網(wǎng)格以獲得可靠的結(jié)果。對(duì)于大型油氣藏礦場(chǎng)尺度的模擬，精細(xì)的網(wǎng)格往往意味著高昂的代價(jià)和冗長的模擬機(jī)時(shí)，其模擬機(jī)時(shí)甚至于超出現(xiàn)有計(jì)算機(jī)的模擬能力。

化學(xué)驅(qū)的驅(qū)替過程往往同時(shí)存在十多個(gè)組分和油、氣、水、乳液等多個(gè)相態(tài)。過程中的各種效應(yīng)高度耦合，存在較強(qiáng)的非線性相互作用。同時(shí)，其驅(qū)替效率又很大程度依賴于局部的化學(xué)劑濃度(如表面活性劑濃度)和乳液相態(tài)變化。故化學(xué)驅(qū)模擬過程對(duì)網(wǎng)格尺寸的敏感性較高，需要精細(xì)的網(wǎng)格以獲得精確的模擬結(jié)果[2-4]。

本文首先綜述用以提高油藏?cái)?shù)值模擬精度的3種主流方法：動(dòng)態(tài)網(wǎng)格加密、粗化技術(shù)和高精度數(shù)值模擬格式。對(duì)于化學(xué)驅(qū)，特別是三元復(fù)合驅(qū)ASP(Alkaline Surfactant Polymer)問題，總結(jié)了前人在動(dòng)態(tài)網(wǎng)格加密和粗化技術(shù)上的成果，并提出了第3種有效的解決方案：使用高精度格式進(jìn)行模擬。其次，利用化學(xué)驅(qū)模擬軟件UTCHEM，對(duì)不同網(wǎng)格尺寸下的一維、二維ASP模型進(jìn)行模擬。發(fā)現(xiàn)模擬結(jié)果對(duì)使用的網(wǎng)格尺寸高度敏感，對(duì)于礦場(chǎng)尺度而言，異常精細(xì)的網(wǎng)格才能獲得相對(duì)收斂的計(jì)算結(jié)果。通過比較模擬結(jié)果，特別是同一時(shí)刻下化學(xué)劑濃度的分布，發(fā)現(xiàn)模似精度對(duì)網(wǎng)格尺寸具有高敏感性是因?yàn)榇蟪叽缇W(wǎng)格對(duì)于化學(xué)組分及化學(xué)劑濃度存在過度的“人工稀釋”效應(yīng)。最后，根據(jù)這一發(fā)現(xiàn)，本文選用高精度格式以更加精確地刻畫網(wǎng)格內(nèi)部組分及化學(xué)劑的濃度分布。數(shù)值模擬測(cè)試表明，高精度格式可以大大提高模擬精度，或者滿足一定精度所需的網(wǎng)格尺寸要求。這對(duì)于降低模擬計(jì)算成本、提高模擬精度有著巨大的幫助。使用高階格式的化學(xué)驅(qū)模擬將為油田礦場(chǎng)尺度大型化學(xué)驅(qū)數(shù)值模擬提供有效的技術(shù)解決方案。

1 油藏?cái)?shù)值模擬精度改進(jìn)的主流方法

主流的商業(yè)或?qū)W術(shù)型油藏?cái)?shù)值模擬器都基于基本的一階迎風(fēng)格式有限差分或有限體積方法[11-12]。其特點(diǎn)為，在網(wǎng)格間流量計(jì)算時(shí)使用上游(迎風(fēng)向)網(wǎng)格內(nèi)平均的壓力、飽和度、組分濃度等物理量來計(jì)算多相達(dá)西定律所決定的流量。一階迎風(fēng)格式的網(wǎng)格內(nèi)物理量完全平均造成較大的模擬誤差。這在復(fù)雜技術(shù)提高采收率過程中尤為顯著，比如火燒油層[9]、混相、非混相氣驅(qū)和凝析氣開采[7-8]等。

改進(jìn)基于一階迎風(fēng)格式的油藏?cái)?shù)值模擬精度有3種主流方法：動(dòng)態(tài)網(wǎng)格加密、粗化技術(shù)和高精度數(shù)值模擬格式。

動(dòng)態(tài)網(wǎng)格加密方法是根據(jù)當(dāng)時(shí)動(dòng)態(tài)的流場(chǎng)信息，僅在油藏中物理量變化劇烈的區(qū)域或重點(diǎn)研究部分使用加密的網(wǎng)格，而其他部位仍使用粗網(wǎng)格的模擬方法。該方法的優(yōu)勢(shì)是能顯著降低計(jì)算模擬成本。缺點(diǎn)是編程實(shí)現(xiàn)較為復(fù)雜，實(shí)際效果高度依賴于使用經(jīng)驗(yàn)，即如何設(shè)定恰當(dāng)?shù)膮?shù)以命令程序在特定區(qū)域內(nèi)加密而在其他區(qū)域內(nèi)粗化。目前，CMG STARS軟件已經(jīng)具備熱采模擬的動(dòng)態(tài)網(wǎng)格加密功能[11]。李建芳[3]等實(shí)驗(yàn)研究了化學(xué)驅(qū)數(shù)值模擬的動(dòng)態(tài)網(wǎng)格加密技術(shù)。

在油藏模擬中利用擬屬性等方法，把細(xì)網(wǎng)格模擬粗化為粗網(wǎng)格模擬,并對(duì)粗網(wǎng)格模擬結(jié)果進(jìn)行矯正的技術(shù)稱為粗化技術(shù)。比如，利用擬相滲曲線方法，可以將黑油問題或多層水驅(qū)過程用粗網(wǎng)格進(jìn)行模擬并保持類似的模擬精度[14-15]。同樣，利用擬分流函數(shù)方法，水氣交替混相氣驅(qū)問題可以使用粗網(wǎng)格進(jìn)行模擬[13]。在火燒油藏模擬中，燃燒反應(yīng)的燃料值作為關(guān)鍵屬性被用于大尺度網(wǎng)格的粗化模擬[9]。擬屬性方法同樣也被運(yùn)用到了化學(xué)驅(qū)的模擬中[10]。同樣地，擬屬性還被用于粗化模型物質(zhì)平衡中網(wǎng)格間流量的計(jì)算。

高階格式也被用于精細(xì)模擬提高采收率過程。如圖1所示，二階格式有限體積方法中，物理量在每個(gè)網(wǎng)格當(dāng)中進(jìn)行線性重構(gòu)(三階格式為二次多項(xiàng)式)。當(dāng)進(jìn)行流量計(jì)算時(shí)，我們代入迎風(fēng)方向的網(wǎng)格邊界處的物理量(如各相飽和度或組分濃度)。目前，求解雙曲型對(duì)流方程(例如油藏模擬中基于達(dá)西定律的每個(gè)組分的物質(zhì)平衡方程)主流的高精度格式為基于Total Variance Diminishing(TVD)概念的格式[5-6]。常見的包括二階TVD格式和三階Essentially Non-Oscillatory(ENO)格式。高階格式對(duì)于混相、非混相氣驅(qū)及凝析氣開采中精確捕捉間斷或激波展現(xiàn)出獨(dú)特的優(yōu)勢(shì)[7-8]。對(duì)于化學(xué)驅(qū)，特別是三元復(fù)合驅(qū)ASP問題，高階精度格式也已在模擬軟件UTCHEM中實(shí)現(xiàn)[12]。本文即在探討該方法對(duì)于降低計(jì)算成本和提高精度的顯著作用。

圖1 一階與二階格式的網(wǎng)格內(nèi)物理量分布：左圖的一階格式網(wǎng)格內(nèi)物理量為平均分布；右圖的二階格式網(wǎng)格內(nèi)物理量為線性分布Fig. 1 Physical property distribution in the fi rst order and second order fi nite difference schemes, with the top fi gure showing fl at distribution of averaged physical properties in fi rst order scheme, and the fi gure below showing piece-wise linear distribution of the physical properties in second order scheme

2 化學(xué)驅(qū)模擬精度的網(wǎng)格敏感性分析

2.1 一維模型模擬精度的網(wǎng)格敏感性分析

建立一維巖芯尺度模型，模型長0.8789 m，寬0.1 m，厚0.1 m，地層原始?jí)毫?.1 MPa，含油飽和度為0.6171。在模型前端置一口注入井，定流量注入，注入過程主要包括4個(gè)階段(詳細(xì)注入液組分及每個(gè)階段的注入速率見表1)： 第1階段注入速度為0.491 mL/min、累計(jì)注入0.679 PV(Pore Volume)；第2階段注入速度為0.052 mL/min、累計(jì)注入0.3 PV的化學(xué)段塞(段塞化學(xué)組成為體積分?jǐn)?shù)2%的表面活性劑、質(zhì)量分?jǐn)?shù)為0.15聚合物、5.11×10-3mol/L的碳酸鈉溶液)；第3階段注入速度為0.052 mL/min、累計(jì)注入1.05 PV；第4階段0.052 mL/min、累計(jì)注入1.7 PV。整個(gè)過程總注入孔隙體積為3.739 PV。模型末端置一口生產(chǎn)井，定井底壓力(0.1 MPa)開采，地層孔隙度為0.1988，x、y、z方向絕對(duì)滲透率為236 mD、236 mD、23.6 mD，地層水黏度0.995 mPa·s，密度1.005 g/cm3，原油黏度24.3 mPa·s，原油密度為0.84 g/cm3，體積系數(shù)為1.101，表面活性劑臨界膠束濃度CMC(Critical Micelle Concentration)為 0.0001 vol%(體積分?jǐn)?shù))。分6類網(wǎng)格尺寸模擬，分別在x方向劃分10、20、40、80、160、240個(gè)網(wǎng)格進(jìn)行先水驅(qū)后ASP三元化學(xué)驅(qū)模擬。在此6種模擬方案中模型的網(wǎng)格數(shù)量不同，其他模擬參數(shù)均一致，采用組分模擬軟件UTCHEM進(jìn)行模擬。

不同網(wǎng)格數(shù)的油藏采收模擬過程見圖2，在不同網(wǎng)格數(shù)模擬下油藏累計(jì)采收率出現(xiàn)差異，最細(xì)網(wǎng)格(240網(wǎng)格)的最終采收率達(dá)到了76.13%，而最粗網(wǎng)格(10網(wǎng)格)的最終采收率僅為63.09%，絕對(duì)采收率相差13.04%。在此模型中，布置240個(gè)網(wǎng)格已經(jīng)達(dá)到足夠精確，如果以240網(wǎng)格的模擬采收率作為基準(zhǔn)，則10網(wǎng)格模擬結(jié)果與240網(wǎng)格模擬結(jié)果相對(duì)誤差達(dá)到17.13%。事實(shí)上在對(duì)礦場(chǎng)尺度油藏模擬時(shí)，即使在此例范圍內(nèi)布置最粗網(wǎng)格(10網(wǎng)格)，由于其占據(jù)巨大的計(jì)算內(nèi)存及計(jì)算量以至于模擬不能正常實(shí)現(xiàn)，這無疑使得模擬結(jié)果更加偏離實(shí)際情況。此外，在圖2中，前期水驅(qū)過程中無論網(wǎng)格尺寸怎樣變化，各油藏模擬采收過程高度一致并與水驅(qū)階段的采收率相同；而在后期ASP三元化學(xué)驅(qū)中，模擬油藏采收率均出現(xiàn)大幅增加，但在不同網(wǎng)格數(shù)模擬采收過程中出現(xiàn)了差異。具體表現(xiàn)為網(wǎng)格尺寸越精細(xì)，油藏最終采收率越高。這說明化學(xué)驅(qū)可以大幅度提高采收率，并且化學(xué)驅(qū)較水驅(qū)對(duì)網(wǎng)格敏感性更高，需要更細(xì)網(wǎng)格進(jìn)行模擬。在此例中，若規(guī)定相對(duì)誤差小于5%，則此模型中采用80網(wǎng)格以上，采收率穩(wěn)定并趨于收斂。

本文進(jìn)一步驗(yàn)證時(shí)間離散誤差對(duì)于此類型三元復(fù)合驅(qū)模擬精度的影響。針對(duì)此模型進(jìn)行時(shí)間步長的敏感性分析。UTCHEM軟件通過設(shè)定模擬中最大和最小CFL數(shù)(Courant–Friedrichs–Lewy)來控制時(shí)間步長。通過改變CFL數(shù)范圍，來設(shè)定不同的時(shí)間步長，以比較不同的時(shí)間離散誤差。這里針對(duì)一維10個(gè)網(wǎng)格的模型進(jìn)行了3組模擬。試驗(yàn)a、b與自動(dòng)調(diào)節(jié)步長實(shí)驗(yàn)c的CFL值范圍分別為：(a)CFL最大值和最小值均為0.0001；(b)CFL最大值和最小值等于0.1；(c)CFL最大值為0.1、最小值為0.0001(在此區(qū)間內(nèi)自動(dòng)調(diào)節(jié))。3個(gè)模型采收過程如圖3所示，a、b、c模型最終油收率分別為64.32%、63.25%、63.09%，差異均小于2%。說明時(shí)間步長對(duì)結(jié)果影響很小，時(shí)間離散誤差對(duì)于此類化學(xué)驅(qū)問題相比空間離散誤差較小，本文不做重點(diǎn)考慮。本文中其他所有一維模型均采用CFL最大值為0.1、最小值為0.0001的自動(dòng)時(shí)間步長設(shè)置。

表1 不同階段(不同累計(jì)注入孔隙體積)下的注入液組分及注入速率Table 1 Composition of injectant and injection speed during different stages, measured by cumulative pore volume injected

為探索化學(xué)驅(qū)較水驅(qū)對(duì)網(wǎng)格具有高敏感度的原因，不同網(wǎng)格數(shù)模擬條件下，在化學(xué)驅(qū)過程中的中間時(shí)刻(驅(qū)替時(shí)間為0.735 d)時(shí)，選取表面活性劑體積分?jǐn)?shù)(圖4)以及剩余油飽和度的分布圖進(jìn)行分析(圖5)。如圖4所示，當(dāng)驅(qū)替時(shí)間為0.735 d時(shí)，與細(xì)網(wǎng)格相比，粗網(wǎng)格表面活性劑濃度出現(xiàn)了“人工稀釋”現(xiàn)象，即粗網(wǎng)格模擬表面活性劑體積分?jǐn)?shù)較細(xì)網(wǎng)格模擬值在峰值處低，但在靠近井口處高。這是粗網(wǎng)格模擬計(jì)算中的數(shù)值彌散現(xiàn)象，這種現(xiàn)象直接影響了化學(xué)驅(qū)驅(qū)替過程中油水界面張力，另一方面可能使得化學(xué)驅(qū)流體相態(tài)發(fā)生變化，導(dǎo)致化學(xué)驅(qū)替過程中驅(qū)油效率降低。驅(qū)替效率可以具體表現(xiàn)為該時(shí)刻下原油剩余飽和度的分布差異如圖5所示。在驅(qū)替前緣經(jīng)過的地層，細(xì)網(wǎng)格剩余原油飽和度較粗網(wǎng)格分布更低，即細(xì)網(wǎng)格驅(qū)油效率高，粗網(wǎng)格驅(qū)油效率低。粗網(wǎng)格出現(xiàn)的化學(xué)劑“人工稀釋”現(xiàn)象導(dǎo)致化學(xué)驅(qū)替過程累計(jì)采收率較低。

圖2 采用不同網(wǎng)格數(shù)的一維ASP模擬的采收過程Fig. 2 1D ASP reservoir recovery process simulation using different number of grid blocks

由此模擬結(jié)果可發(fā)現(xiàn)，在同一模型下，網(wǎng)格數(shù)量影響了化學(xué)驅(qū)模擬的精度，網(wǎng)格數(shù)量越多或網(wǎng)格劃分越細(xì)其化學(xué)劑精度越高，模擬結(jié)果更精確。

圖3 采用不同的時(shí)間步長(即不同的CFL值)下的一維ASP驅(qū)數(shù)值模擬采收率(空間離散為10個(gè)網(wǎng)格)Fig. 3 The reservoir recovery of one-dimensional ASP fl ood using different time step sizes (different CFL values)

2.2 二維模型模擬精度的網(wǎng)格敏感性分析

圖4 采用不同網(wǎng)格數(shù)進(jìn)行一維ASP模擬得到的表面活性劑體積分?jǐn)?shù)分布(0.735 d)Fig. 4 The volume fraction distribution of surfactant in simulation under different number of grid blocks in one dimensional ASP fl ood (0.735 days)

建立二維油藏模型，長1 m，寬1 m，厚0.1 m，巖石性質(zhì)、流體性質(zhì)及注入方式與2.1中一維模型參數(shù)保持一致，注入井與生產(chǎn)井呈對(duì)角分布，上、下對(duì)角分別布置1口注入井、1口生產(chǎn)井，采用4種不同網(wǎng)格尺寸對(duì)同一模型進(jìn)行模擬，即在長、寬、高方 向 分 別 布 置 5×5×1、10×10×1、20×20×1、30×30×1網(wǎng)格數(shù)。油藏采收率模擬結(jié)果如圖6，結(jié)果證明二維模擬與一維模擬具有一致性，即化學(xué)驅(qū)模擬精度高度依賴網(wǎng)格尺寸大小或網(wǎng)格數(shù)目，網(wǎng)格尺寸越小，網(wǎng)格越細(xì)，化學(xué)驅(qū)模擬結(jié)果越精確。

圖5 采用不同網(wǎng)格數(shù)進(jìn)行模擬得到一維ASP模擬的原油飽和度分布(0.735 d)Fig. 5 Oil saturation distribution in simulation under different number of grid blocks in one dimensional ASP fl ood (0.735 days)

圖6 采用不同網(wǎng)格數(shù)下的二維ASP模擬的采收過程Fig. 6 Simulated recovery factor to pore volume injected in two dimensional ASP model by using different number of grid blocks

與一維模型分析類似，取二維模型化學(xué)驅(qū)驅(qū)替中間時(shí)刻7.35 d時(shí)，表面活性劑體積分?jǐn)?shù)及剩余油飽和度分布如圖7、圖8所示，隨著網(wǎng)格數(shù)的減少，網(wǎng)格內(nèi)部的表面活性劑體積分?jǐn)?shù)呈現(xiàn)出與一維模型一致的現(xiàn)象，即表面活性劑呈現(xiàn)“人工稀釋”現(xiàn)象。出現(xiàn)這種“人工稀釋”現(xiàn)象時(shí)，此時(shí)的剩余油飽和度分布(圖8)顯示，粗網(wǎng)格模擬在驅(qū)替前緣經(jīng)過的地層的剩余油飽和度，要高于細(xì)網(wǎng)格模擬。這說明在化學(xué)驅(qū)驅(qū)替階段的低網(wǎng)格數(shù)目模擬的驅(qū)油效率降低了。這種現(xiàn)象是造成粗網(wǎng)格模擬油藏累計(jì)采收率降低的原因。

圖7 使用5×5×1、10×10×1、20×20×1網(wǎng)格時(shí)的二維ASP模擬的表面活性劑體積分?jǐn)?shù)分布(7.35 d)Fig. 7 Surfactant volume fraction in 2D ASP process at 7.35 days using 5×5×1, 10×10×1, 20×20×1 grid system

圖8 使用5×5×1、10×10×1、20×20×1網(wǎng)格時(shí)的二維ASP模擬的原油飽和度分布(7.35 d)Fig. 8 Oil saturation distribution in 2D ASP process at 7.35 days using 5×5×1, 10×10×1, 20×20×1 grid system

與一維時(shí)間離散誤差分析類似，在此針對(duì)二維5×5×1模型采用3種不同CFL值設(shè)定以檢驗(yàn)時(shí)間步長對(duì)于二維模型的影響。3種模型CFL數(shù)分別為：(a)CFL最大值和最小值均為0.0001；(b)CFL最大值和最小值均為0.1；(c)CFL最大值為0.1、最小值為0.0001(在此區(qū)間內(nèi)自動(dòng)調(diào)節(jié))。油藏采收情況如圖9所示，a、b、c模型最終采收率分別為52.87%、52.81%、52.89%，三者采收率差異小于1%。由此得知，二維模型結(jié)果對(duì)于時(shí)間步長不敏感。本文中其他所有二維模型均采用CFL最大值為0.1、最小值為0.0001的自動(dòng)時(shí)間步長設(shè)置。

圖9 在不同CFL值(即時(shí)間步長)下的二維ASP模型的采收率(5×5×1網(wǎng)格)Fig. 9 Recovery of two-dimensional 5×5×1 ASP model under different CFL values

3 高階格式化學(xué)驅(qū)模擬

3.1 化學(xué)驅(qū)高階格式求解方法

首先，考慮一維多組分多相多孔介質(zhì)流動(dòng)的控制方程：

這里Ci和Fi分別是第i個(gè)組分的濃度和流量。uj為第j相的達(dá)西流動(dòng)速度：

當(dāng)采用有限體積方法求解該問題時(shí)，第k個(gè)網(wǎng)格的離散化的控制方程為：

Δt和Δx分別為時(shí)間和空間步長，k?1 2和k+1 2分別為第i個(gè)組分在第k個(gè)網(wǎng)格左側(cè)和右側(cè)網(wǎng)格邊界的流動(dòng)通量。當(dāng)采用一階迎風(fēng)格式時(shí)，網(wǎng)格內(nèi)的物理量采用常值進(jìn)行重構(gòu)(零階多項(xiàng)式)。一階格式，當(dāng)?shù)趉個(gè)網(wǎng)格右側(cè)網(wǎng)格邊界的流量向右及向左時(shí)，k+1 2處的通量分別為時(shí)：

當(dāng)采用二階格式時(shí)，網(wǎng)格內(nèi)的物理量采用線性重構(gòu)(一階多項(xiàng)式)。當(dāng)?shù)趉個(gè)網(wǎng)格右側(cè)網(wǎng)格邊界的流量向右時(shí)，其處的通量為：

二階格式往往伴隨著間斷處的數(shù)值震蕩，需要加入一定的數(shù)值黏性以消除震蕩，基于Total Variance Diminishing概念的二階TVD格式隨之而生[5-6]。其核心理念為設(shè)計(jì)一系列數(shù)值格式，以控制物理量的Total Variance不隨時(shí)間而放大：

常見的二階TVD格式采取以下形式：

φ(rk)為Flux Limiter函數(shù)。當(dāng)φ(rk)=1時(shí)，格式退化為普通二階格式。這里的兩種Flux Limiter函數(shù)分別為常見的Min mod函數(shù)和Fromm函數(shù)。

高階格式(如三階的ENO格式)的構(gòu)建與二階格式類似，即在考慮迎風(fēng)的情況下對(duì)網(wǎng)格內(nèi)部的物理量進(jìn)行二階及更高階的多項(xiàng)式重構(gòu)，并計(jì)算網(wǎng)格邊界處k?1 2、k+1 2處的通量。同樣加入修正項(xiàng)以保證對(duì)間斷的正確捕捉和防止震蕩產(chǎn)生。高階格式帶松弛的Total Variance限制條件為(n為格式的階數(shù))：三維空間的多組分多相油藏?cái)?shù)值模擬與之類似?？紤]迎風(fēng)的情況下，一階格式、二階格式和三階格式在網(wǎng)格內(nèi)分別對(duì)物理量進(jìn)行空間3個(gè)方向上的常值、線性和二次多項(xiàng)式重構(gòu)。利用類似的修正方法，高階格式在間斷附近的震蕩可被消除，最終達(dá)到精確描述流動(dòng)過程的目的。

3.2 高階格式一維化學(xué)驅(qū)模擬

分別用二階、三階格式對(duì)章節(jié)2.1中一維巖心模型進(jìn)行模擬運(yùn)算，高階格式運(yùn)算使用UTCHEM軟件內(nèi)置的功能，其油藏采收過程如圖10所示。

結(jié)果表明：二階、三階格式模擬采收率結(jié)果比一階格式更接近細(xì)網(wǎng)格模擬結(jié)果。具體表現(xiàn)為：采用20×1×1網(wǎng)格的三階格式的模擬結(jié)果與240×1×1一階格式的結(jié)果十分接近，也可以說在采收率精度上兩者相當(dāng)；三階格式模擬采收率則較二階格式更接近細(xì)網(wǎng)格模擬結(jié)果。對(duì)化學(xué)驅(qū)中間時(shí)刻(驅(qū)替時(shí)間為0.735 d)的表面活性劑體積分?jǐn)?shù)分布(圖11)以及殘余油飽和度分布(圖12)進(jìn)行分析，發(fā)現(xiàn)出現(xiàn)了同樣的現(xiàn)象。如20網(wǎng)格三階格式下表面活劑和剩余油飽和度分布與240網(wǎng)格一階格式模擬結(jié)果基本一致。此模型中，在達(dá)到240網(wǎng)格一階格式結(jié)果精度的條件下，三階格式模擬僅需使用20網(wǎng)格。

圖10 使用不同格式及不同網(wǎng)格數(shù)下一維ASP模型的采收過程Fig. 10 One-dimensional ASP recovery process by using different number of grid blocks and different numerical schemes

圖11 使用不同格式及不同網(wǎng)格數(shù)下的一維ASP模擬的表面活性劑分布(0.735 d)Fig. 11 The volume fraction distribution of surfactant in one-dimensional ASP model by using different numerical schemes and different number of grid blocks (0.735 days)

圖12 使用不同格式及不同網(wǎng)格數(shù)下的一維ASP模擬的原油飽和度分布(0.735 d)Fig. 12 Oil saturation distribution in one-dimensional ASP model by using different numerical schemes and different number of grid blocks (0.735 days)

3.3 高階格式二維化學(xué)驅(qū)模擬

對(duì)二維模型使用不同階格式進(jìn)行模擬運(yùn)算，其模型參數(shù)與章節(jié)2.2中二維模型一致，模擬結(jié)果油藏采收率見圖13，從圖13中可以看出二維油藏模擬中，高階格式模擬運(yùn)算情況下，其油藏采收率精度更接近于細(xì)網(wǎng)格精度，如10×10×1三階格式運(yùn)算油藏采收率為60.59%，這與細(xì)網(wǎng)格30×30×1一階格式運(yùn)算油藏采收率60.53%幾乎相等，這證實(shí)在二維模型中高階格式一樣可以改善化學(xué)驅(qū)精度，實(shí)現(xiàn)低網(wǎng)格數(shù)目高精度的目的。

圖13 使用不同格式及不同網(wǎng)格數(shù)下的二維ASP模擬采收過程Fig. 13 Simulated two-dimensional ASP recovery process using different difference schemes and different number of grid blocks

取驅(qū)替時(shí)間為7.35天時(shí)，不同網(wǎng)格尺寸配合不同階格式模擬的表面活性劑體積分?jǐn)?shù)分布情況如圖14所示，高階格式表面活性劑精度要高于一階格式模擬精度，網(wǎng)格越細(xì)，表面活性劑精度越高。存在高階粗網(wǎng)格精度與細(xì)網(wǎng)格精度相當(dāng)?shù)呐R界粗網(wǎng)格數(shù)，如10×10×1三階格式粗網(wǎng)格表面活性劑精度可以代表30×30×1一階格式細(xì)網(wǎng)格表面活性劑精度。

3.4 高階格式三維礦場(chǎng)尺度化學(xué)驅(qū)模擬

建立礦場(chǎng)尺度多井組化學(xué)驅(qū)油藏?cái)?shù)值模擬模型。模擬區(qū)塊長200 m，寬200 m，厚度10 m。模型縱向分4層，第1層在x、y、z方向的滲透率分別為236 mD、236 mD、23.6 mD，第2到第4層的滲透率分別為第1層滲透率的0.8、0.6、0.4倍。地層壓力為20.7 MPa。該區(qū)塊采用4個(gè)5點(diǎn)井網(wǎng)進(jìn)行開發(fā)，包括4口注入井和9口生產(chǎn)井。注入井采用定流量控制，分三個(gè)階段。第一階段為水驅(qū)，注入量為141.4 m3/d、累計(jì)注入0.679 PV。第二階段為注入化學(xué)段塞(段塞含有體積分?jǐn)?shù)為0.02的表面活性劑、質(zhì)量分?jǐn)?shù)為0.15的聚合物及濃度為5.11×10-3mol/L的碳酸鈉)，注入量為14.9 m3/d、累計(jì)注入0.3 PV。第三階段為(注入液主要成分包括濃度為2.5×10-3mol/L氯化鈉及濃度為0.04×10-3mol/L氯化鎂)，注入量為14.9 m3/d、累計(jì)注入1.05 PV。以上3個(gè)階段的總注入液量為2.0368 PV。生產(chǎn)井則采用井底壓力控制生產(chǎn)(6.89 MPa)。該區(qū)塊的地層流體屬性與巖石物性與章節(jié)2.1中的案例一致。分別使用11×11×4、21×21×4、41×41×4三組由粗到細(xì)的網(wǎng)格來模擬三元復(fù)合驅(qū)過程。各井的位置分布如圖15(以41×41×4網(wǎng)格為例)所示。對(duì)11×11×4、21×21×4網(wǎng)格模型分別采用一階、二階、三階高階格式進(jìn)行數(shù)值模擬。而細(xì)網(wǎng)格模型41×41×4則采用一階格式，模擬采用化學(xué)驅(qū)模擬器UTCHEM。

該礦場(chǎng)尺度化學(xué)驅(qū)模型采用不同網(wǎng)格及不同格式進(jìn)行模擬得到的采收率結(jié)果如圖16所示，與一維、二維的模擬類似，當(dāng)采用不同格式進(jìn)行礦場(chǎng)尺度模型模擬時(shí)，粗網(wǎng)格配合高階格式的運(yùn)算精度可與細(xì)網(wǎng)格配合低階格式的運(yùn)算精度相當(dāng)，例如21×21×4網(wǎng)格配合三階格式的采收率為56.45%與41×41×4網(wǎng)格配合一階格式的采收率為56.51%一致。

為探究采收過程差異原因，選用模擬時(shí)間為143 d時(shí)，不同格式下油藏的第1層的表面活性劑體積分?jǐn)?shù)分布如圖17所示。原油飽和度的分布如圖18所示。21×21×4網(wǎng)格配合三階格式與41×41×4網(wǎng)格在采收率、表面活性劑分布及原油飽和度分布精度一致。即高階格式下粗網(wǎng)格模擬的表面活性劑體積分?jǐn)?shù)分布與低階格式細(xì)網(wǎng)格相當(dāng)。其他屬性如原油飽和度分布也一致。這證明了礦場(chǎng)尺度下，高階格式可以有效提高粗網(wǎng)格的模擬精度。

圖14 使用不同格式及不同網(wǎng)格數(shù)下二維ASP模擬的表面活性劑體積分?jǐn)?shù)的分布(7.35 d)Fig. 14 Surfactant volume fraction distribution in two-dimensional ASP flooding under different difference schemes and different number of grid blocks (7.35 days)

圖15 礦場(chǎng)尺度多井組三元復(fù)合驅(qū)模型中4口注入井和9口生產(chǎn)井的井位分布及41×41×4網(wǎng)格模型中驅(qū)替時(shí)刻為343 d時(shí)的原油飽和度分布(上箭頭代表生產(chǎn)井，下箭頭代表注入井)Fig. 15 Well locations including four injection wells and nine production wells in multi-well group fi led-scale ASP model and distribution of oil saturation at the time of 343 days in 41×41×4 grid model

圖17 使用不同格式及不同網(wǎng)格數(shù)下三維礦場(chǎng)尺度模型的表面活性劑體積分?jǐn)?shù)分布(143 d，圖中P代表生產(chǎn)井，I代表注入井)Fig. 17 The volume fraction distribution of surfactant in three-dimensional fi led-scale model using different difference schemes and different grid numbers (143 days, P represents a production well and I represents a injection well)

圖18 使用不同格式及不同網(wǎng)格數(shù)下三維礦場(chǎng)尺度模型的原油飽和度分布圖(143 d，圖中P代表生產(chǎn)井，I代表注入井)Fig. 18 The distribution of oil saturation using different difference schemes and different grid numbers for 3D fi led-scale model(143 days, P represents a production well and I represents a injection well)

4 計(jì)算成本對(duì)比

在此對(duì)使用不同網(wǎng)格尺寸及不同格式的2.1中一維模型模擬運(yùn)算CPU運(yùn)行時(shí)間進(jìn)行了統(tǒng)計(jì)和對(duì)比(圖19)。章節(jié)3.2已證明一維ASP模型使用20個(gè)網(wǎng)格配合三階格式與使用240個(gè)網(wǎng)格配合一階格式的運(yùn)算精度基本類似。從圖19中發(fā)現(xiàn)一維ASP模擬使用20個(gè)網(wǎng)格和三階格式的CPU占用時(shí)間為5.32秒，而同等精度的使用240個(gè)網(wǎng)格和一階格式的CPU占用時(shí)間為269.8秒。所以在本例當(dāng)中，在運(yùn)算精度相同的條件下，運(yùn)用粗網(wǎng)格三階格式模擬占用的CPU時(shí)間僅為細(xì)網(wǎng)格一階格式CPU時(shí)間的2%。前者極大縮短了模型的運(yùn)算時(shí)間、降低了計(jì)算成本，并且粗網(wǎng)格模擬由于網(wǎng)格數(shù)低，占用計(jì)算機(jī)內(nèi)存小，同樣減輕了計(jì)算機(jī)的硬件配置要求。這都提高了化學(xué)驅(qū)油藏?cái)?shù)值模擬效率。

在二維模型中，高階格式也能得到相同的效果(圖20)。使用10×10×1網(wǎng)格配合三階格式占用的CPU時(shí)間為46.5秒，而等同精度的30×30×1網(wǎng)格配合一階格式所占用CPU時(shí)間則為2676秒。本例中，細(xì)網(wǎng)格一階格式模擬所占用的CPU時(shí)間為粗網(wǎng)格三階格式的57.5倍，網(wǎng)格數(shù)是粗網(wǎng)格模擬的9倍。高階格式具有非常明顯的優(yōu)勢(shì)。

三維礦場(chǎng)尺度模型由于其區(qū)塊尺寸大，網(wǎng)格數(shù)量更多，往往需要更高的模擬運(yùn)行時(shí)間。在本文中三維礦場(chǎng)尺度多井組化學(xué)驅(qū)模型被離散為11×11×4、21×21×4、41×41×4三組不同的網(wǎng)格，并采用不同的格式進(jìn)行模擬。但CPU運(yùn)行時(shí)間有十分大的差異，如圖21所示。章節(jié)3.4中已證明在三維礦場(chǎng)尺度模型油藏?cái)?shù)值模擬中，21×21×4網(wǎng)格配合三階格式與41×41×4配合一階格式具有等同精度，此例中粗網(wǎng)格(21×21×4)配合高階格式CPU占用時(shí)間為885.6秒，而細(xì)網(wǎng)格(41×41×4)配合一階格式CPU占用時(shí)間高達(dá)7264.3秒，在同等模擬精度條件下，粗網(wǎng)格配合高階格式數(shù)值模擬占用計(jì)算機(jī)CPU運(yùn)行時(shí)間更少，計(jì)算機(jī)運(yùn)行成本更低，在高精度化學(xué)驅(qū)油藏?cái)?shù)值模擬中顯現(xiàn)出了極大的優(yōu)勢(shì)。

圖19 使用不同網(wǎng)格數(shù)配合不同格式下的一維ASP模型模擬占用的CPU時(shí)間對(duì)比，此處分別標(biāo)注了同一模型同等模擬精度，但網(wǎng)格數(shù)和格式不同的模擬算例，及其相應(yīng)的計(jì)算成本的差別Fig. 19 The CPU run time comparison for one-dimensional model using different grid numbers and different differential schemes, with simulations having the same model, and using different grid numbers and numerical schemes, but of similar accuracies marked, which corresponds to different CPU runtime

圖20 使用不同網(wǎng)格數(shù)及不同格式下的二維ASP模型模擬占用的CPU時(shí)間對(duì)比，此處分別標(biāo)注了同一模型同等模擬精度，但網(wǎng)格數(shù)和格式不同的模擬算例，及其相應(yīng)的計(jì)算成本的差別Fig. 20 The CPU run time comparison for two-dimensional models using different grid numbers and different differential schemes, with simulations having the same model, and using different grid numbers and numerical schemes, but of similar accuracies marked, which corresponds to different CPU runtime

圖21 使用不同網(wǎng)格數(shù)和不同格式下的三維礦場(chǎng)尺度模型模擬占用的CPU時(shí)間對(duì)比，此處分別標(biāo)注了同一模型同等模擬精度，但網(wǎng)格數(shù)和格式不同的模擬算例，及其相應(yīng)的計(jì)算成本的差別Fig. 21 The Comparison of CPU time simulated by 3D fi led-scale model with different grid numbers and difference schemes,with simulations having the same model, and using different grid numbers and numerical schemes, but of similar accuracies marked, which corresponds to different CPU runtime

5 結(jié)論

(1)通過不同網(wǎng)格尺寸下一維、二維化學(xué)驅(qū)ASP模擬，我們發(fā)現(xiàn)ASP三元復(fù)合驅(qū)模擬精度對(duì)于網(wǎng)格尺寸具有高度敏感，且在同樣網(wǎng)格尺寸下ASP模擬的精度遠(yuǎn)低于水驅(qū)模擬精度。

(2)化學(xué)驅(qū)模擬中，當(dāng)使用大尺寸網(wǎng)格時(shí)，組分和化學(xué)劑濃度(特別是表面活性劑濃度)出現(xiàn)過度的“人工稀釋”現(xiàn)象。這一數(shù)值假象造成了模擬中精度的降低，即粗網(wǎng)格情況下驅(qū)油效率的下降。

(3)采用基于TVD原理的高階(二階TVD、三階ENO)格式進(jìn)行離散后的化學(xué)驅(qū)數(shù)值模擬，其精度有很大的提升。這為大型礦場(chǎng)尺度化學(xué)驅(qū)油藏?cái)?shù)值模擬提供了有效的提高精度并且降低計(jì)算成本的方法。

(4)在一維、二維、三維模型中分別測(cè)試了不同網(wǎng)格尺寸配合不同格式模型的CPU運(yùn)行時(shí)間，無論哪種模型都充分顯示了粗網(wǎng)格配合高階格式較細(xì)網(wǎng)格配合低階格式所占用的計(jì)算機(jī)CPU時(shí)間更少，計(jì)算精度更高。

[1] 韓大匡. 油藏?cái)?shù)值模擬基礎(chǔ)[M]. 北京: 石油工業(yè)出版社, 1993: 2-9. [HAN D K. Fundamentals of numerical reservoir simulation[M].Beijing: Petroleum Industry Press, 1993: 2-9.]

[2] 劉皖露, 馬德勝, 王強(qiáng), 等. 化學(xué)驅(qū)數(shù)值模擬技術(shù)[J]. 大慶石油學(xué)院學(xué)報(bào), 2012, 36(3): 72-78. [LIU W L, MA D S, WANG Q, et al.Numerical simulation for chemical fl ooding[J]. Journal of Daqing Petroleum Institute, 2012, 36 (3): 72-78.]

[3] 李建芳, 袁士義, 宋杰, 等. 化學(xué)驅(qū)驅(qū)替前緣動(dòng)態(tài)追蹤數(shù)值模擬研究[J]. 石油勘探與開發(fā), 2004, 31(B11): 55-58. [LI J F, YUAN S Y, SONG J, et al. Reservoir simulation of chemical fl ooding by dynamically tracking displacement fronts[J]. Petroleum Exploration and Development, 2004, 31 (B11): 55-58.]

[4] 袁士義. 注化學(xué)劑驅(qū)油數(shù)值模擬 (應(yīng)用部分)[J]. 石油學(xué)報(bào), 1989, 10(3): 68-76. [YUAN S Y. Note Chemical agent fl ooding numerical simulation (application part)[J]. Acta Petrolei Sinica, 1989, 10 (3): 68-76.]

[5] SWEBY P K. High resolution schemes using fl ux limiters for hyperbolic conservation laws[J]. SIAM journal on numerical analysis,1984, 21(5): 995-1011.

[6] HARTEN A. High resolution schemes for hyperbolic conservation laws[J]. Journal of computational physics, 1983, 49(3): 357-393.

[7] MALLISON B T, GERRITSEN M G, JESSEN K, et al. High order upwind schemes for two-phase, multicomponent flow[J]. SPE Journal, 2005, 10(03): 297-311.

[8] JESSEN K, GERRITSEN M G, MALLISON B T. High-resolution prediction of enhanced condensate recovery processes[J]. SPE Journal, 2008, 13(02): 257-266.

[9] NISSEN A, ZHU Z, KOVSCEK A, et al. Upscaling kinetics for fi eld-scale in-situ-combustion simulation[J]. SPE Reservoir Evaluation& Engineering, 2015, 18(02): 158-170.

[10] KOYASSAN VEEDU F, DELSHAD M, POPE G A. Scaleup methodology for chemical fl ooding[C]//SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, 2010.

[11] CHRISTENSEN J R, DARCHE G, DECHELETTE B, et al. Applications of dynamic gridding to thermal simulations[C]//SPE international thermal operations and heavy oil symposium and western regional meeting. Society of Petroleum Engineers, 2004.

[12] DELSHAD M, ASAKAWA K, POPE G A, et al. Simulations of chemical and microbial enhanced oil recovery methods[C]//SPE/DOE Improved Oil Recovery Symposium. Society of Petroleum Engineers, 2002.

[13] CHRISTIE M A, MANSFIELD M, KING P R, et al. A renormalisation-based upscaling technique for WAG fl oods in heterogeneous reservoirs[C]//SPE Reservoir Simulation Symposium. Society of Petroleum Engineers, 1995.

[14] HEARN C L. Simulation of strati fi ed water fl ooding by pseudo relative permeability curves[J]. Journal of Petroleum Technology, 1971,23(07): 805-813.

[15] STONE H L. Rigorous Black Oil Pseudo Functions[C]// SPE Symposium on Reservoir Simulation. Society of Petroleum Engineers,1991.

[16] CHEN Z, HUAN G, MA Y. Computational methods for multiphase fl ows in porous media[C]. Computational Science and Engineering Series, Vol. 2, SIAM, Philadelphia, 2006.

[17] CHEN Z, MA Y, CHEN G. A sequential numerical chemical compositional simulator[J]. Transport in Porous Media, 2007, 68: 389-411.