賈英，湯昌福，黃磊

（1.中國石化石油勘探開發(fā)研究院，北京 100083；2.安徽省煤田地質(zhì)局勘查研究院，安徽 合肥 230088；3.中國石油勘探開發(fā)研究院，北京 100083）

基于迎風(fēng)差分格式的TVD聚合物驅(qū)數(shù)學(xué)模擬

賈英1，湯昌福2，黃磊3

（1.中國石化石油勘探開發(fā)研究院，北京 100083；2.安徽省煤田地質(zhì)局勘查研究院，安徽 合肥 230088；3.中國石油勘探開發(fā)研究院，北京 100083）

聚合物驅(qū)替過程中存在飽和度和聚合物質(zhì)量濃度2個(gè)間斷，在2個(gè)間斷之間有“油墻”形成。文中分別運(yùn)用特征線法（MOC）、一階迎風(fēng)差分（SPU）方法和顯式全變差遞減差分（TVD）方法對聚合物驅(qū)替過程進(jìn)行求解，通過與特征線解進(jìn)行對比，考察了SPU方法和TVD方法的數(shù)值模擬效果。研究表明：SPU方法的精度低，數(shù)值耗散嚴(yán)重，對2個(gè)間斷特別是聚合物質(zhì)量濃度間斷的捕捉能力較差；高階精度TVD方法分辨力高，但在“油墻”處會(huì)產(chǎn)生非物理“振蕩”。利用迎風(fēng)的思想，對傳統(tǒng)TVD方法進(jìn)行了改進(jìn)。數(shù)值實(shí)驗(yàn)表明，改進(jìn)的TVD方法在保持高精確度的同時(shí)，還能有效地防止非物理“振蕩”的產(chǎn)生，具有較好的穩(wěn)定性。

聚合物驅(qū)；特征線法；一階迎風(fēng)格式；TVD；非物理“振蕩”

0 引言

聚合物驅(qū)能夠有效改變驅(qū)替流度比，擴(kuò)大掃油面積，提高波及系數(shù)和采出程度。近年來，聚合物驅(qū)提高采收率技術(shù)得到了大力推廣和應(yīng)用［1-9］。

聚合物驅(qū)是一個(gè)復(fù)雜的物理過程，在特定的情形下，可以運(yùn)用特征線法MOC（Method of Characteristics）進(jìn)行解析求解［2-3］，然而，在大部分情況下，聚合物驅(qū)都需要借助于數(shù)值方法進(jìn)行研究［3，7］。在油藏?cái)?shù)值模擬中，對流占優(yōu)滲流方程主要采用單點(diǎn)上游迎風(fēng) SPU（Single Point Upstream）方法求解［10］。然而，SPU方法只是一階精度，具有較強(qiáng)數(shù)值耗散作用，驅(qū)替過程出現(xiàn)的激波（間斷）可能會(huì)被抹平［11-16］。高階精度差分方法雖可提高數(shù)值模擬的精度，但直接使用容易在間斷附近產(chǎn)生非物理“振蕩”［14］?？傋儾钸f減（TVD）方法通過通量限制因子控制模擬過程中總變差的變化，在獲取高階精度的同時(shí)，還能有效防止非物理“振蕩”的出現(xiàn)［13-14］；因而，TVD方法在油藏?cái)?shù)值模擬中已得到推廣和應(yīng)用［11-15］。

本文首先運(yùn)用MOC方法求解一維聚合物驅(qū)過程，通過與MOC解進(jìn)行對比，進(jìn)一步分析了聚合物驅(qū)SPU和TVD數(shù)值解的效果，提出了改進(jìn)的TVD數(shù)值通量的重構(gòu)方法，在保證了傳統(tǒng)TVD方法的高階精度的同時(shí)，還防止了非物理“振蕩”的出現(xiàn)。

1 基本數(shù)學(xué)模型

聚合物驅(qū)主要通過增加水相黏度減小驅(qū)替流度比來提高采收率。在聚合物驅(qū)替過程中，還伴有吸附、水相滲透率下降等現(xiàn)象［1-3］。本文主要考慮了水相黏度的變化及巖石吸附對聚合物驅(qū)替的影響。并且假設(shè)：均質(zhì)不可壓等溫流動(dòng)；油水兩相流動(dòng)，聚合物溶于水中；驅(qū)替過程對流占優(yōu)，不考慮擴(kuò)散作用，且忽略重力和毛細(xì)管力對驅(qū)替過程的影響。聚合物溶液的黏度是聚合物質(zhì)量濃度的函數(shù)，可用修正的水相黏度方程表示：

聚合物驅(qū)的吸附為等溫吸附過程，可以用線性吸附或Langmuir吸附表示。為研究方便，假設(shè)線性吸附過程，即聚合物吸附是聚合物質(zhì)量濃度的線性函數(shù)：

一維聚合物驅(qū)數(shù)學(xué)模型的水相連續(xù)方程和聚合物質(zhì)量濃度方程分別為

將式（3）、式（4）無量綱化，得：

將式（6）、式（7）二元雙曲型方程組寫成矩陣形式：

式（8）的邊界條件和初始條件分別為

聚合物溶液的黏度隨聚合物質(zhì)量濃度呈線性變化。關(guān)于聚合物質(zhì)量濃度的多項(xiàng)式H（C）為

研究中，取a1=1 cm3/mg，a2=9 cm3/mg，μw=1 mPa·s，μo=10 mPa·s。

采用Corey型相滲曲線，即：

取Swi=0.2，Sor=0.3，Kw=0.6，Ko=1.0，no=2，nw=2。

設(shè)注入聚合物溶液質(zhì)量濃度Cin=1 mg/cm3，吸附相關(guān)常數(shù)D=0.1。油、水相對滲透率曲線及水相分流量曲線見圖1。

圖1 油、水相對滲透率及水相分流量曲線

式（8）和式（9）組成了一維聚合物驅(qū)的Riemann問題，可以采用MOC求解。

2 一維聚合物驅(qū)特征線法

Riemann問題的解可以由簡單波描述，簡單波由稀疏波和激波（或接觸間斷）組成。式（8）存在2個(gè)特征值，分別為

圖2（τ=0.1）為一維聚合物驅(qū)Riemann問題MOC模擬結(jié)果?？梢钥闯?，聚合物驅(qū)過程中存在兩族波，分別為S族和C族。其中：S族波為飽和度間斷，聚合物質(zhì)量濃度C為常數(shù)；C族波由接觸間斷和稀疏波組成，在接觸間斷處，聚合物質(zhì)量濃度C和飽和度Sw都是間斷的，但其特征速度λc為常數(shù)。飽和度間斷和聚合物質(zhì)量濃度間斷分別決定了聚合物驅(qū)見水和見聚的時(shí)間。在兩族波之間存在“油墻”，飽和度保持為常數(shù)。

圖2 MOC模擬飽和度和質(zhì)量濃度空間分布情況

一維聚合物驅(qū)Riemann問題的MOC解具體可參考文獻(xiàn)［2-3］。該計(jì)算結(jié)果與巖心實(shí)驗(yàn)CT掃描結(jié)果相一致（見圖3）。

圖3 水飽和度曲線及CT分析數(shù)據(jù)［16］

3 一維聚合物驅(qū)有限差分法

當(dāng)初始條件和注入條件保持穩(wěn)定時(shí)，一維聚合物驅(qū)可以運(yùn)用特征線法（MOC）快速準(zhǔn)確地求解。然而，實(shí)際應(yīng)用過程中，飽和度和聚合物質(zhì)量濃度的初始分布一般都不均勻，且注入條件也經(jīng)常發(fā)生變化，此時(shí)常需要采用數(shù)值方法求解，所采用的數(shù)值方法應(yīng)能夠準(zhǔn)確地捕捉聚合物驅(qū)中出現(xiàn)的物理現(xiàn)象。

式（8）的顯式差分方法守恒差分格式可以表示為

不同的通量限制因子（Φ（r））構(gòu)成不同的差分方法。

3.1 一階迎風(fēng)（SPU）差分方法

在油藏?cái)?shù)值模擬中，由于具有明確的物理特點(diǎn)，經(jīng)常采用單點(diǎn)上游加權(quán)迎風(fēng)格式SPU［10］求解，即數(shù)值流通量采用其上游的網(wǎng)格值計(jì)算：

圖4（τ=0.1）為不同網(wǎng)格數(shù)時(shí)SPU解與MOC解析解之間的對比?？梢钥闯觯壕W(wǎng)格數(shù)較少時(shí)，聚合物驅(qū)過程中的2個(gè)間斷——特別是質(zhì)量濃度間斷——都被不同程度抹平了，SPU的數(shù)值模擬結(jié)果嚴(yán)重偏離精確解；隨著網(wǎng)格數(shù)的增加，SPU數(shù)值解趨于MOC解，此時(shí)為獲得一個(gè)較為精確的差分解，需要更多網(wǎng)格。

圖4 一維聚合物驅(qū)不同網(wǎng)格數(shù)SPU差分模擬結(jié)果

上述結(jié)果表明，SPU方法對聚合物驅(qū)過程中的2個(gè)間斷特別是質(zhì)量濃度間斷的捕捉能力較弱。這是因?yàn)椋航佑|間斷與激波不同，它不具有“自銳利”（Selfsharpening）行為，精確解的一些特征被差分方法具有的數(shù)值耗散掩蓋，需要采用其他方法（例如粒子追蹤方法、高階精度差分方法等［11，17］）才能改進(jìn)模擬的結(jié)果。

3.2 高分辨率顯式全變差遞減（TVD）差分方法

TVD差分方法能提高數(shù)值模擬的的精度，然而直接采用該方法容易在間斷附近產(chǎn)生非物理“振蕩”。為了定量研究“振蕩”情形，定義時(shí)間步n的總變差［11］為

TVD差分方法通過限制總變差的增加，防止非物理“振蕩”的產(chǎn)生，即：

此時(shí)式（14）中的通量限制因子Φ為梯度r的函數(shù)，使得在模擬區(qū)域的不同部分需采用不同的差分格式。目前，學(xué)者們提出了各種限制因子，如Min mod限制因子、Fromm限制因子、Van Leer限制因子以及Leonard限制因子等［8］。本文采用Van Leer限制因子的TVD差分方法研究一維聚合物驅(qū)過程：

圖5（τ=0.1）為一維聚合物驅(qū)TVD差分方法模擬結(jié)果。通過與SPU解以及MOC解進(jìn)行對比可以發(fā)現(xiàn)，聚合物驅(qū)數(shù)值模擬的效果得到顯著改善，在較少的網(wǎng)格數(shù)時(shí)就能獲得較為精確的數(shù)值解 （N=400）。因此，TVD差分方法捕捉間斷的能力要強(qiáng)于SPU方法，聚合物驅(qū)差分方法應(yīng)優(yōu)先選擇TVD差分方法，保證在較少的網(wǎng)格時(shí)得到更加精確的數(shù)值解。

圖5 一維聚合物驅(qū)TVD差分模擬結(jié)果

然而，從圖5還可以發(fā)現(xiàn)，TVD差分方法在提高精度的同時(shí)，卻在 “油墻”處產(chǎn)生了非物理 “振蕩”。Thiele等研究發(fā)現(xiàn)，增加時(shí)間差分精度、控制CFL的大小或者采用耗散性較強(qiáng)的限制因子（如Min mod）可以減緩“振蕩”，但并不能達(dá)到消除“振蕩”的目的。

3.3 迎風(fēng)TVD差分方法

傳統(tǒng)TVD差分方法在“油墻”處會(huì)出現(xiàn)非物理“振蕩”?；谟筒?cái)?shù)值模擬中常用的迎風(fēng)思想，本文提出了新的數(shù)值流通量的計(jì)算方法，以抑制“振蕩”的發(fā)生。在新的TVD差分方法中，飽和度方程中的數(shù)值流通量采用迎風(fēng)格式，而質(zhì)量濃度方程中的總變差通過質(zhì)量濃度的TVD來限制，即：

聚合物質(zhì)量濃度方程中數(shù)值通量依式（21）計(jì)算：

圖6（τ=0.1）給出了改進(jìn)的TVD解（U-TVD）與MOC解及傳統(tǒng)TVD解的對比?？梢园l(fā)現(xiàn)，迎風(fēng)TVD差分方法在保持TVD差分方法高階精度的同時(shí)，還有效地抑制了“油墻”處非物理“振蕩”的出現(xiàn)。

4 結(jié)論

1）一階迎風(fēng)（SPU）方法模擬聚合物驅(qū)能力較弱。網(wǎng)格數(shù)較少時(shí)，數(shù)值解嚴(yán)重偏離精確解；雖然可通過增加網(wǎng)格數(shù)改善數(shù)值模擬結(jié)果，但增加的計(jì)算代價(jià)巨大。

2）TVD差分方法能以較少網(wǎng)格獲得較精確的數(shù)值解，聚合物驅(qū)數(shù)值模擬應(yīng)盡量采用高分辨率差分方法。

3）傳統(tǒng)TVD差分方法會(huì)在飽和度和質(zhì)量濃度間斷之間存在的“油墻”處產(chǎn)生不合理的非物理“振蕩”。基于迎風(fēng)思想，改進(jìn)了TVD數(shù)值通量計(jì)算方法，在保證TVD方法的高精度同時(shí)，還能防止“振蕩”產(chǎn)生，使得TVD方法更為穩(wěn)定。該方法為下一步開展三維聚合物驅(qū)數(shù)值模擬奠定了基礎(chǔ)。

圖6 一維聚合物驅(qū)迎風(fēng)TVD差分模擬結(jié)果

5 符號(hào)注釋

μp，μw，μo分別為聚合物溶液、水相和油相黏度，mPa·s；C為聚合物在水相中的質(zhì)量濃度，mg/cm3；Ca為聚合物吸附質(zhì)量濃度，mg/cm3；Cin為聚合物溶液質(zhì)量濃度，mg/cm3；a為吸附常數(shù)；x為空間坐標(biāo)，m；t為時(shí)間坐標(biāo)，s；ρr，ρw分別為巖石和水相密度，kg/m3；q為注入流量，m3/s；Sw為水飽和度；Swi，Soi分別為束縛水飽和度和殘余油飽和度；φ為孔隙度；A為截面積，m2；f為水相分流量函數(shù) （即水飽和度和聚合物質(zhì)量濃度的函數(shù)）；Krw，Kro分別為水相和油相的相對滲透率；Kw，Ko分別為殘余油飽和度端點(diǎn)水相滲透率和束縛水飽和度端點(diǎn)油相滲透率；τ為注入孔隙體積倍數(shù)；Vp為孔隙體積，m3；L為模型長度，m；xD為無量綱長度；D是吸附相關(guān)常數(shù)；a1，a2為常數(shù) （通過實(shí)驗(yàn)取值）；nw，no分別為水相、油相相對滲透率曲線Coery系數(shù)；n和i分別表示時(shí)間和空間節(jié)點(diǎn)（或網(wǎng)格位置）；N為空間網(wǎng)格數(shù)；和分別為網(wǎng)格i的左、右數(shù)值流通量；Φ（）r為通量限制因子；r為梯度；下標(biāo)分別表示上一網(wǎng)格和下一網(wǎng)格；表示在網(wǎng)格位置的質(zhì)量濃度差比值。

［1］NEEDHAM R B，DOE P H.Polymer flooding review［J］.Journal of Petroleum Technology，1987，39（12）：1503-1507.

［2］POPE GARY A.The application of fractional flow theory to enhanced oil recovery［J］.SPE Journal 1980，20（20）：191-205.

［3］PATTON J T，COATS K H，COLEGROVE G T.Prediction of polymer flood performance［R］.SPE 2546，1971.

［4］蔣明.齊40塊稠油油藏聚合物驅(qū)可行性室內(nèi)實(shí)驗(yàn)研究［J］.斷塊油氣田，1998，5（1）：45-49.

［5］宋考平，楊二龍，王錦梅，等.聚合物驅(qū)提高驅(qū)油效率機(jī)理及驅(qū)油效果分析［J］.石油學(xué)報(bào)，2004，25（3）：71-74.

［6］方文超.陸上稠油油田多輪次聚合物驅(qū)提高采收率方法［J］.斷塊油氣田，2015，22（5）：619-622.

［7］蘆文生.綏中36-1油田聚合物驅(qū)數(shù)值模擬研究［J］.中國海上油氣，2002，16（5）：333-340.

［8］周叢叢.大慶油田一、二類油層聚合物驅(qū)注采指標(biāo)變化規(guī)律［J］.斷塊油氣田，2015，22（5）：610-613.

［9］侯健，李振泉，王玉斗，等.考慮擴(kuò)散和吸附作用的聚合物驅(qū)替過程滲流數(shù)值模擬［J］.計(jì)算物理，2003，20（3）：239-244.

［10］MALLISON B T，GERRITSEN M G，JESSEN K，et al.Higher-order upwind schemes for two-phase，multi-component Flow［R］.SPE 79691，2005.

［11］THIELE M R，EDWARDS M.Physically based higher order godunov schemes for compositional simulation［R］.SPE 66403，2001.

［12］LEER V B.Towards the Ultimate conservative difference Scheme，V-A second-order sequel to Godunov′s method［J］.Journal of Computational Physics，1979，32（1）：101-136.

［13］HARTEN A.High resolution schemes for hyperbolic conservation laws［J］.Journal of Computational Physics，1983，49（3）：357-393.

［14］RUBIN B，EDWARDS M G.Extension of the TVD mid-point scheme to higher order accuracy in time［R］.SPE 25265，1993.

［15］BRENIER Y，JAFFRé J.Upstream differencing for multiphase flow in reservoir simulation［J］.Siam Journal on Numerical Analysis，1991，28（3）：685-696.

［16］SU S.Laboratory investigation of polymer-flooding，ASP-flooding，and foam-flooding in heavy oil system［D］.Calgary：University of Calgary，2012.

［17］JUANES R，LIE K A.Numerical modeling of multiphase first-contact miscible flows.Ⅱ.Front-tracking/streamline simulation［J］.Transport Porous Media，2008，72（1）：97-120.

（編輯 李宗華）

Mathematical simulation of polymer flooding with upstream-TVD solution

JIA Ying1,TANG Changfu2,HUANG Lei3
(1.Petroleum Exploration and Production Research Institute,SINOPEC,Beijing 100083,China;2.Exploration Research Institute, Anhui Provincial Bureau of Coal Geology,Anhui 230088,China;3.Research Institute of Petroleum Exploration&Development, PetroChina,Beijing 100083,China)

For polymer flooding,there are two discontinuities of saturation and concentration,between which"oil bank"exists.The polymer flooding was investigated with method of characteristics(MOC),single point upstream(SPU)finite difference method and explicit total variation diminishing(TVD)finite difference method respectively.Based on the solution of MOC,performance of simulation results of both SPU and TVD were evaluated.It is suggested TVD exhibits a better numerical performance while the SPU scheme suffers numerical diffusion badly.The conventional TVD method is modified on the calculation of numerical flux,so that non-physical oscillation can be prevented on the"oil bank"between the two discontinuities.Numerical tests demonstrate that the upstream-TVD can achieve better stability than the conventional TVD method.

polymer flooding;MOC;SPU;TVD;non-physical oscillation

TE357.46+1

A

10.6056/dkyqt201703017

2016-12-12；改回日期：2017-03-09。

賈英，女，1980年生，高級(jí)工程師，博士，油氣田開發(fā)工程專業(yè)。E-mail：jiaying.syky@sinopec.com。

賈英，湯昌福，黃磊.基于迎風(fēng)差分格式的TVD聚合物驅(qū)數(shù)學(xué)模擬［J］.斷塊油氣田，2017，24（3）：373-377.

JIA Ying，TANG Changfu，HUANG Lei.Mathematical simulation of polymer flooding with upstream-TVD solution［J］.Fault-Block Oil&Gas Field，2017，24（3）：373-377.