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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3.1 移動(dòng)最小二乘法

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

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

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

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

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

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

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

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

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

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

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

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

其中，

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

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

2 算例分析

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

圖2 一注四采模型與節(jié)點(diǎn)布置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]

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

圖3 營養(yǎng)物和微生物混合溶液不同注入量下代謝產(chǎn)物質(zhì)量濃度對(duì)比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ǔ)之上，對(duì)影響無網(wǎng)格法計(jì)算精度和計(jì)算效率的布點(diǎn)方案、節(jié)點(diǎn)影響域和權(quán)函數(shù)分別進(jìn)行參數(shù)敏感性分析，進(jìn)而驗(yàn)證該方法在微生物驅(qū)數(shù)值模擬求解中的正確性。

3.1 不同布點(diǎn)密度分析

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

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

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

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

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

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

不同的權(quán)函數(shù)能使節(jié)點(diǎn)在其影響域里所產(chǎn)生的影響隨權(quán)函數(shù)不同而改變，其對(duì)最終結(jié)果的精確度也會(huì)產(chǎn)生一定影響。本文對(duì)微生物驅(qū)模型選用了不同的權(quán)函數(shù)，選用9×9 的布點(diǎn)方案，影響域中的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ù)計(jì)算含水率變化曲線Fig.7.Variations of water cut with different weight functions in the microbial flooding model

4 結(jié)論

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

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

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

符號(hào)注釋

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

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

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

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

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

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

dj——xj點(diǎn)到x點(diǎn)的距離，m；

dmj——影響域半徑，m；

Dm——微生物在水相中的擴(kuò)散系數(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é)點(diǎn)；

Sw——含水飽和度；

t——時(shí)間，d；

Δt——時(shí)間步長，d；

u——局部最佳近似；

vt——流速，m/s；

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

α——罰函數(shù)；

φ——孔隙度；

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