基于無網(wǎng)格法的微生物驅(qū)數(shù)值模擬

曹琳，修建龍，黃立信，黃峰，徐云峰，趙輝，盛廣龍

（1.中國科學(xué)院大學(xué)，北京 100049；2.中國科學(xué)院滲流流體力學(xué)研究所，河北廊坊 065007；3.中國石油勘探開發(fā)研究院，北京 100083；4.中國石油新疆油田分公司重油開發(fā)公司，新疆克拉瑪依 834000；5.長江大學(xué)石油工程學(xué)院，武漢 430100）

經(jīng)歷了傳統(tǒng)的一次采油和二次采油之后，油藏中仍然有大約2/3 的原油剩余[1-3]。現(xiàn)有提高采收率技術(shù)中，微生物驅(qū)提高采收率技術(shù)是一項施工簡單、費用低廉、安全環(huán)保、不傷害油層、返排液無需處理、潛力巨大的三次采油技術(shù)[4]。微生物驅(qū)是利用微生物及其代謝產(chǎn)物作用于油層，降低原油黏度和油水界面張力，從而提高油藏中原油的流動性，最終提高采收率[5-6]。目前，微生物驅(qū)機理研究主要采用室內(nèi)實驗結(jié)合數(shù)值模擬的方法。室內(nèi)實驗方面，文獻［7］進行了微生物驅(qū)室內(nèi)實驗，引入了微生物因子對殘余油飽和度的影響。文獻［8］通過實驗研究數(shù)據(jù)，擬合出微生物質(zhì)量濃度與原油黏度之間的關(guān)系，并以此建立了微生物降黏模型。室內(nèi)實驗?zāi)M可以較為真實地反映地下微生物驅(qū)機理，然而由于實驗耗時長、成本高，很難應(yīng)用到實際生產(chǎn)井模擬分析。為此，學(xué)者們普遍采用數(shù)值模擬進行微生物驅(qū)模擬，文獻［9］基于菌體流動及吸附實驗，發(fā)現(xiàn)采油過程中，微生物在多孔介質(zhì)運移中篩分作用比吸附作用更加突出，進而完善了微生物運移數(shù)學(xué)模型；文獻［10］分析了微生物質(zhì)量濃度、營養(yǎng)液質(zhì)量濃度、段塞結(jié)構(gòu)、油層滲透率等因素對驅(qū)油效率的影響，并采用全隱式有限差分法對其所建立數(shù)學(xué)模型進行求解。微生物驅(qū)機理的研究已經(jīng)相對完善，但是在微生物驅(qū)數(shù)學(xué)模型的求解上，主要采用的是有限差分法，此類方法以網(wǎng)格劃分為前提，在計算過程中經(jīng)常需要進行網(wǎng)格重構(gòu)，勢必會對計算速度產(chǎn)生較大影響，同時此方法前后處理也比較麻煩。

1977 年，文獻［11］在解決無邊界天體物理問題時提出了無網(wǎng)格法，其本質(zhì)是利用一組散布在問題域及域邊界上的節(jié)點，構(gòu)造場函數(shù)近似表達式的一種計算方法。文獻［12］用伽遼金無網(wǎng)格法對單相穩(wěn)定滲流問題進行了研究，研究結(jié)果表明，無網(wǎng)格法無需劃分網(wǎng)格，能夠消除劃分網(wǎng)格所帶來的困難，解決許多基于網(wǎng)格方法不能解決的問題。文獻［13］針對油水兩相油藏滲流問題，分析伽遼金無網(wǎng)格法和加權(quán)最小二乘無網(wǎng)格法的基本原理，并將無網(wǎng)格法引入到數(shù)值試井中，建立復(fù)雜斷塊油藏油水兩相無網(wǎng)格試井解釋模型。大量基于無網(wǎng)格法的研究表明，無網(wǎng)格法在計算中完全脫離了網(wǎng)格，前處理工作相對于有限差分法簡單了許多[14-16]。加權(quán)最小二乘無網(wǎng)格法直接使用最小二乘法建立系統(tǒng)的變分原理，消除了控制方程在邊界節(jié)點處的殘差，無需進行GAUSS積分，計算量小。

由于微生物驅(qū)機理復(fù)雜，涉及組分相對較多，使用有限差分法等數(shù)值模擬方法求解時存在計算不穩(wěn)定、計算效率低的缺點[17]。本文將無網(wǎng)格法拓展到微生物驅(qū)數(shù)值模擬中，通過考慮微生物生長/死亡、營養(yǎng)消耗、產(chǎn)物生成、化學(xué)趨向性、對流擴散、吸附、解吸等物理化學(xué)過程，建立了微生物驅(qū)提高采收率數(shù)學(xué)模型，采用加權(quán)最小二乘無網(wǎng)格法進行了求解，分析無網(wǎng)格法在微生物驅(qū)數(shù)值模擬中的應(yīng)用效果。

1 數(shù)學(xué)模型的建立與求解

1.1 微生物驅(qū)物理模型

本文通過構(gòu)建一注一采的一維油藏模型進行微生物驅(qū)機理研究（圖1）。圖中W1 井為注水井，P1 井為采油井，采用均勻布點，序號為節(jié)點數(shù)。

圖1 微生物驅(qū)物理模型Fig.1.Physical model of microbial flooding

模型基本假設(shè)條件：①油藏均質(zhì)且儲集層內(nèi)流體為油、水兩相；②油、水是微可壓縮流體，流體混合無體積變化；③油藏等溫；④微生物及其營養(yǎng)物和代謝產(chǎn)物均只溶于水相中。

1.2 微生物驅(qū)數(shù)學(xué)模型

為定量描述微生物驅(qū)這一過程，建立的微生物驅(qū)數(shù)學(xué)模型不僅需要考慮滲流場中油和水的運移，而且還需包含微生物及其營養(yǎng)物和代謝產(chǎn)物的運移。滲流場中的油、水流動方程滿足達西定律和質(zhì)量守恒定律，其加權(quán)最小二乘無網(wǎng)格法的計算形式可參考文獻［18］。

考慮微生物及其營養(yǎng)物和代謝產(chǎn)物的對流、擴散、吸附、解吸、沉淀等作用，結(jié)合微生物反應(yīng)動力學(xué)，微生物及其營養(yǎng)物和代謝產(chǎn)物的運移方程[19]為

假定井以定產(chǎn)量進行生產(chǎn)，設(shè)定初始含水飽和度，該模型能夠通過設(shè)置不同條件來模擬外源和內(nèi)源微生物驅(qū)，其區(qū)別在于油藏原始微生物質(zhì)量濃度是否為0。對于邊界條件，生物場的質(zhì)量濃度邊界條件均設(shè)為注入常數(shù)。隨著微生物的生長和繁殖，油藏滲透率、原油黏度和相對滲透率會隨之發(fā)生改變，同時這些物性參數(shù)的變化也反映出微生物驅(qū)機理。

1.3 無網(wǎng)格法求解

加權(quán)殘量法可作為無網(wǎng)格法的理論基礎(chǔ)，本文應(yīng)用移動最小二乘法來構(gòu)造近似函數(shù)，然后利用最小二乘法對控制方程進行離散，建立對應(yīng)的無網(wǎng)格法，即加權(quán)最小二乘無網(wǎng)格法。

1.3.1 移動最小二乘法

設(shè)在求解域Ω內(nèi)有場函數(shù)u（x）以及一組隨機分布的離散節(jié)點xj（j=1，2，…，n），用Ωj表示節(jié)點xj的影響域，在二維問題中，影響域常取作圓形或矩形。使用移動最小二乘法進行全局近似，求解域內(nèi)的任意節(jié)點x有近似函數(shù)：

權(quán)函數(shù)是移動最小二乘法近似中的重要組成部分，沒有理論上的具體規(guī)則，帶有任意性，一般滿足權(quán)函數(shù)ω（r）非負，r=dj/dmj，dj=︱x-xj︱，dmj是xj的影響域半徑，本文中影響域半徑取dmj=scale×s［1］，s［1］為節(jié)點xj與距其最近的節(jié)點距離，scale為大于等于1的乘子。

最小二乘法能夠用來求解具有任何類型的偏微分方程，從變分原理和加權(quán)殘量法上講這類方法是來建立泛函的，且其權(quán)函數(shù)或稱為檢驗函數(shù)就是殘量本身。首先，令ΠΩ=再取最小值，即對邊界條件運用更為簡單、有效的本質(zhì)邊界條件處理方案即罰函數(shù)法。罰函數(shù)的功能，是對非可行點或企圖穿越邊界而逃離可行域的點賦予一個極大值，即將有約束最優(yōu)化問題轉(zhuǎn)化為無約束最優(yōu)化問題。

對具體問題采用上述基于移動最小二乘法的近似函數(shù)，并采用上述積分格式的離散點求和形式便得到了加權(quán)最小二乘無網(wǎng)格法[20]。

1.3.2 加權(quán)最小二乘無網(wǎng)格法求解微生物驅(qū)

設(shè)求解域Ω及邊界Γ內(nèi)有n個離散節(jié)點，營養(yǎng)物質(zhì)量濃度近似函數(shù)通過移動最小二乘法來進行全局近似，如下式：

其中，N=[N1(x)N2(x)…Nn(x)]為形函數(shù)；pN=[CN1CN2…CNn]T為營養(yǎng)物質(zhì)量濃度分布的矩陣形式。

將（3）式代入（1）式中，可得營養(yǎng)物運移方程的殘量方程：

因此營養(yǎng)物運移方程的加權(quán)平方和公式為

將（5）式取最小值，即：

上式通常需要寫成離散形式，即：

其中，n1為第一類邊界條件Γ上的節(jié)點總數(shù)，δRN為殘量方程RN的變分形式，寫成矩陣形式如下：

其中，

（10）式求出下一時刻微生物質(zhì)量濃度后，代謝產(chǎn)物質(zhì)量濃度變化速率如下式：

代謝產(chǎn)物運移方程的加權(quán)最小二乘無網(wǎng)格法求解方法與求解營養(yǎng)物和微生物類似，這里不再累述。

2 算例分析

基于上文建立的加權(quán)最小二乘無網(wǎng)格法求解微生物驅(qū)數(shù)學(xué)模型，并拓展到二維問題的求解上，設(shè)計了如下算例，算例井網(wǎng)布置采用一注四采模型（圖2）。W1為注水井，P1井、P2井、P3井和P4井為采油井。正方形油藏邊長為100 m，層厚為40 m，滲透率為100 mD，孔隙度為20%。把方形油藏用9×9 個節(jié)點進行離散，應(yīng)用加權(quán)最小二乘無網(wǎng)格法進行求解，節(jié)點影響域中的scale 為3，選用GAUSS 權(quán)函數(shù)和二次基函數(shù)，時間步長取Δt=1 d，生產(chǎn)制度注采均為47.67 m3/d。本文模型與文獻［19］建立的模型中驅(qū)油機理大致相同，物性參數(shù)如表1所示。

本文算例采用如下生產(chǎn)制度：在注入井一次注水1.65 PV 至含水率95%，后續(xù)以連續(xù)注入的方式注入0.30 PV 營養(yǎng)物和微生物的混合溶液（營養(yǎng)物質(zhì)量濃度為2.495 mg/mL，微生物質(zhì)量濃度為1.275 mg/mL），最終水驅(qū)至含水率98%，微生物菌體在多孔介質(zhì)中持續(xù)吸收營養(yǎng)物，生長、繁殖并代謝出相關(guān)產(chǎn)物，隨著注入的持續(xù)進行，代謝產(chǎn)物的分布會發(fā)生變化。

圖2 一注四采模型與節(jié)點布置Fig.2.Model with one injection well and four producing wells,and deployment of mesh points

表1 加權(quán)最小二乘無網(wǎng)格法求解微生物驅(qū)數(shù)學(xué)模型物性參數(shù)[19-21]Table 1.Physical parameters of the mathematical model for microbial flooding solved by the weighted least square meshless method[19-21]

為了驗證本文方法的正確性與計算效率，將上述算例計算結(jié)果與有限差分法進行比較。有限差分法計算所用參數(shù)和本文所用參數(shù)保持一致，采用12.5 m×12.5 m×40.0 m 的網(wǎng)格作為單個計算單元。兩者計算所得到的代謝產(chǎn)物質(zhì)量濃度分布如圖3 所示，含水率變化與采收率變化曲線如圖4 所示，加權(quán)最小二乘無網(wǎng)格法與有限差分法的計算結(jié)果基本吻合，說明了本方法的正確性與有效性。在加權(quán)最小二乘無網(wǎng)格法求解的節(jié)點單元數(shù)比有限差分法求解的要多15 個的情況下，計算效率反而提高了62%。在總注入量為1.66 PV 時，油藏中營養(yǎng)物和微生物混合溶液比較少，產(chǎn)生的代謝產(chǎn)物也少，主要集中在注入井附近。當總注入量為1.95 PV 時，營養(yǎng)物和微生物混合溶液停止注入，進入后續(xù)水驅(qū)階段后，油藏中的微生物菌體因得不到相應(yīng)的營養(yǎng)物，也將停止產(chǎn)生代謝產(chǎn)物。隨著后續(xù)水驅(qū)階段的進行，代謝產(chǎn)物質(zhì)量濃度峰值也逐步向采出端移動，且代謝產(chǎn)物質(zhì)量濃度峰值也因吸附作用而逐漸下降。

圖3 營養(yǎng)物和微生物混合溶液不同注入量下代謝產(chǎn)物質(zhì)量濃度對比Fig.3.Mass concentration of metabolites at different injected volumes of the mixed solution of nutrient and microorganism

圖4 微生物驅(qū)和水驅(qū)的含水率變化曲線（a）和采收率變化曲線（b）Fig.4.(a)Water cut curves and(b)recovery curves of microbial flooding and water flooding

3 無網(wǎng)格法敏感性分析

在上述理論和求解方法的基礎(chǔ)之上，對影響無網(wǎng)格法計算精度和計算效率的布點方案、節(jié)點影響域和權(quán)函數(shù)分別進行參數(shù)敏感性分析，進而驗證該方法在微生物驅(qū)數(shù)值模擬求解中的正確性。

3.1 不同布點密度分析

應(yīng)用加權(quán)最小二乘無網(wǎng)格法求解時，選擇合適的布點方案非常重要。本文采用等距離均勻分布，在正方形油藏區(qū)域100 m×100 m，選用了6×6、9×9、21×21和51×51 的布點方案，節(jié)點影響域中的scale 為3，權(quán)函數(shù)公式為GAUSS 函數(shù)，進行微生物驅(qū)數(shù)學(xué)模型的計算，并將計算結(jié)果和有限差分法計算結(jié)果作對比。51×51 布點方案的計算精度與有限差分法最接近，但是當布點數(shù)下降時，誤差并沒有大幅增加（圖5）。但布點密的方案所需的計算時間明顯高于布點稀疏的方案，如51×51布點方案的計算時間大約是9×9布點方案的10倍，分別是53.0 s和5.3 s。而9×9 布點方案和6×6布點方案的計算時間卻基本相同。由此可見，求解微生物驅(qū)模型時，在沒有大幅降低計算精度的情況下，每100 m×100 m的油藏區(qū)域采用9×9 的布點密度能大幅度提升計算效率。

3.2 不同的節(jié)點影響域分析

圖5 微生物驅(qū)不同布點方案含水率變化曲線Fig.5.Water cut curves of different meshing schemes for microbial flooding

加權(quán)最小二乘無網(wǎng)格法自身具有較高的精度，但其結(jié)果還依賴于節(jié)點影響域的選取。其參數(shù)最優(yōu)值和所求解的具體問題有關(guān)，因此本文選取節(jié)點影響域中的scale 為2、3、4 和5，對微生物驅(qū)數(shù)學(xué)模型求解進行了驗證分析。從圖6 中能夠看出隨著節(jié)點影響域的減小，計算精度呈現(xiàn)先增加后降低的趨勢，在scale為3 時，計算精度最高。對計算效率分析發(fā)現(xiàn)，節(jié)點影響域?qū)τ嬎阈什]有太大影響，運行時間為5.5 s左右，而有限差分法運行時間為13.8 s。因此本文選用scale為3來對微生物驅(qū)數(shù)學(xué)模型進行求解，以確保求解的正確性。

圖6 微生物驅(qū)不同節(jié)點影響域含水率分布Fig.6.Water cut distribution in the influence domain of different points

3.3 不同的權(quán)函數(shù)分析

不同的權(quán)函數(shù)能使節(jié)點在其影響域里所產(chǎn)生的影響隨權(quán)函數(shù)不同而改變，其對最終結(jié)果的精確度也會產(chǎn)生一定影響。本文對微生物驅(qū)模型選用了不同的權(quán)函數(shù)，選用9×9 的布點方案，影響域中的scale為3，權(quán)函數(shù)分別用SPLIN、GAUSS、CUBIC 和CSRBF。從圖7中可以看出，采用GAUSS權(quán)函數(shù)的最終結(jié)果精確度最高，說明這種權(quán)函數(shù)分布最適合于微生物驅(qū)數(shù)學(xué)模型的求解。

圖7 微生物驅(qū)不同權(quán)函數(shù)計算含水率變化曲線Fig.7.Variations of water cut with different weight functions in the microbial flooding model

4 結(jié)論

（1）將無網(wǎng)格法引入微生物驅(qū)數(shù)學(xué)模型求解，考慮對流、擴散、微生物繁殖和死亡、生物趨化性、營養(yǎng)物消耗、代謝產(chǎn)物生成、吸附等特性，建立了微生物提高采收率數(shù)學(xué)模型，采用無網(wǎng)格法對滲流場和生物場進行了求解。

（2）將無網(wǎng)格法與有限差分法求解微生物驅(qū)數(shù)學(xué)模型進行了對比驗證，結(jié)果表明，在確保計算正確性的前提下，無網(wǎng)格法計算效率明顯優(yōu)于有限差分法。

（3）對影響無網(wǎng)格法計算精度和計算效率的參數(shù)進行了敏感性分析，得出在邊長為100 m 方形油藏區(qū)域采用9×9 均勻布點方案，影響域為3 倍節(jié)點間距為半徑的圓形和GAUSS 權(quán)函數(shù)時，更適用于微生物驅(qū)數(shù)學(xué)模型的求解，能在確保計算正確性的前提下，提高計算的效率。

符號注釋

aN、bN——營養(yǎng)物的吸附常數(shù)；

A——菌體維持生命時代謝產(chǎn)物的生成速率，1/d；

B——組分代謝產(chǎn)物得率，g/g；

Bw——地層體積系數(shù)；

Ci、Cis——分別為水相和吸附相中各組分的質(zhì)量濃度，g/L；

CN、CNs——分別為營養(yǎng)物質(zhì)量濃度和第s個點的營養(yǎng)物質(zhì)量濃度，g/L；

dj——xj點到x點的距離，m；

dmj——影響域半徑，m；

Dm——微生物在水相中的擴散系數(shù)，m2/d；

i——N（營養(yǎng)物）、m（微生物）和d（代謝產(chǎn)物）；

kc——微生物的化學(xué)趨向系數(shù)，m2/d；

N——無網(wǎng)格法中的形函數(shù)；

qw——日注入量（采出量），m3；

R2——代謝產(chǎn)物質(zhì)量濃度變化速率，g/d；

Rd——組分脫附速率，g/L；

Ri——組分變化速率，g/L；

RN、Rm——分別為營養(yǎng)物和微生物控制方程的殘量方程；

Rr——組分吸附速率，g/L；

RΩ——任意控制方程在求解域內(nèi)的產(chǎn)量方程；

s——無網(wǎng)格法中代表節(jié)點；

Sw——含水飽和度；

t——時間，d；

Δt——時間步長，d；

u——局部最佳近似；

vt——流速，m/s；

Vb——源匯處巖石孔隙體積，m3；

α——罰函數(shù)；

φ——孔隙度；

?！刂品匠痰倪吔鐥l件。