基于高階格式的高精度化學驅(qū)模擬

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

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

0 引言

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

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

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

1 油藏數(shù)值模擬精度改進的主流方法

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

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

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

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

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

圖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 化學驅(qū)模擬精度的網(wǎng)格敏感性分析

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

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

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

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

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

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

圖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ù)量影響了化學驅(qū)模擬的精度，網(wǎng)格數(shù)量越多或網(wǎng)格劃分越細其化學劑精度越高，模擬結(jié)果更精確。

圖3 采用不同的時間步長(即不同的CFL值)下的一維ASP驅(qū)數(shù)值模擬采收率(空間離散為10個網(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ù)進行一維ASP模擬得到的表面活性劑體積分數(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)井呈對角分布，上、下對角分別布置1口注入井、1口生產(chǎn)井，采用4種不同網(wǎng)格尺寸對同一模型進行模擬，即在長、寬、高方 向 分 別 布 置 5×5×1、10×10×1、20×20×1、30×30×1網(wǎng)格數(shù)。油藏采收率模擬結(jié)果如圖6，結(jié)果證明二維模擬與一維模擬具有一致性，即化學驅(qū)模擬精度高度依賴網(wǎng)格尺寸大小或網(wǎng)格數(shù)目，網(wǎng)格尺寸越小，網(wǎng)格越細，化學驅(qū)模擬結(jié)果越精確。

圖5 采用不同網(wǎng)格數(shù)進行模擬得到一維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

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

圖7 使用5×5×1、10×10×1、20×20×1網(wǎng)格時的二維ASP模擬的表面活性劑體積分數(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)格時的二維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

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

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

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

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

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

這里Ci和Fi分別是第i個組分的濃度和流量。uj為第j相的達西流動速度：

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

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

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

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

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

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

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

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

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

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

對二維模型使用不同階格式進行模擬運算，其模型參數(shù)與章節(jié)2.2中二維模型一致，模擬結(jié)果油藏采收率見圖13，從圖13中可以看出二維油藏模擬中，高階格式模擬運算情況下，其油藏采收率精度更接近于細網(wǎng)格精度，如10×10×1三階格式運算油藏采收率為60.59%，這與細網(wǎng)格30×30×1一階格式運算油藏采收率60.53%幾乎相等，這證實在二維模型中高階格式一樣可以改善化學驅(qū)精度，實現(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ū)替時間為7.35天時，不同網(wǎng)格尺寸配合不同階格式模擬的表面活性劑體積分數(shù)分布情況如圖14所示，高階格式表面活性劑精度要高于一階格式模擬精度，網(wǎng)格越細，表面活性劑精度越高。存在高階粗網(wǎng)格精度與細網(wǎng)格精度相當?shù)呐R界粗網(wǎng)格數(shù)，如10×10×1三階格式粗網(wǎng)格表面活性劑精度可以代表30×30×1一階格式細網(wǎng)格表面活性劑精度。

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

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

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

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

圖14 使用不同格式及不同網(wǎng)格數(shù)下二維ASP模擬的表面活性劑體積分數(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 礦場尺度多井組三元復合驅(qū)模型中4口注入井和9口生產(chǎn)井的井位分布及41×41×4網(wǎng)格模型中驅(qū)替時刻為343 d時的原油飽和度分布(上箭頭代表生產(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ù)下三維礦場尺度模型的表面活性劑體積分數(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ù)下三維礦場尺度模型的原油飽和度分布圖(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 計算成本對比

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

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

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

圖19 使用不同網(wǎng)格數(shù)配合不同格式下的一維ASP模型模擬占用的CPU時間對比，此處分別標注了同一模型同等模擬精度，但網(wǎng)格數(shù)和格式不同的模擬算例，及其相應的計算成本的差別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時間對比，此處分別標注了同一模型同等模擬精度，但網(wǎng)格數(shù)和格式不同的模擬算例，及其相應的計算成本的差別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ù)和不同格式下的三維礦場尺度模型模擬占用的CPU時間對比，此處分別標注了同一模型同等模擬精度，但網(wǎng)格數(shù)和格式不同的模擬算例，及其相應的計算成本的差別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)格尺寸下一維、二維化學驅(qū)ASP模擬，我們發(fā)現(xiàn)ASP三元復合驅(qū)模擬精度對于網(wǎng)格尺寸具有高度敏感，且在同樣網(wǎng)格尺寸下ASP模擬的精度遠低于水驅(qū)模擬精度。

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

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

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

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

[2] 劉皖露, 馬德勝, 王強, 等. 化學驅(qū)數(shù)值模擬技術(shù)[J]. 大慶石油學院學報, 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] 李建芳, 袁士義, 宋杰, 等. 化學驅(qū)驅(qū)替前緣動態(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] 袁士義. 注化學劑驅(qū)油數(shù)值模擬 (應用部分)[J]. 石油學報, 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.