裂縫性碳酸鹽巖酸化模擬新方法

羅 源楊乾龍劉平禮趙立強(qiáng)杜 勛徐 彬（.西南石油大學(xué)油氣藏地質(zhì)及開(kāi)發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川成都 60500；.中國(guó)石油長(zhǎng)慶油田分公司第十采油廠，甘肅慶陽(yáng) 74500）

引用格式：羅源，楊乾龍，劉平禮，等.裂縫性碳酸鹽巖酸化模擬新方法［J］.石油鉆采工藝，2015，37（5）：73-77.

裂縫性碳酸鹽巖酸化模擬新方法

羅 源1楊乾龍2劉平禮1趙立強(qiáng)1杜 勛1徐 彬1
（1.西南石油大學(xué)油氣藏地質(zhì)及開(kāi)發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川成都 610500；2.中國(guó)石油長(zhǎng)慶油田分公司第十采油廠，甘肅慶陽(yáng) 745100）

引用格式：羅源，楊乾龍，劉平禮，等.裂縫性碳酸鹽巖酸化模擬新方法［J］.石油鉆采工藝，2015，37（5）：73-77.

摘要：在以往裂縫性碳酸鹽巖酸化模擬的研究中，一般假設(shè)裂縫與基質(zhì)均勻切割且將求解域離散為規(guī)則網(wǎng)格的形式，處理復(fù)雜形態(tài)的裂縫時(shí)會(huì)有較大困難。為了更加真實(shí)地模擬裂縫性碳酸鹽巖的酸液流動(dòng)過(guò)程，采用有限元網(wǎng)格對(duì)其求解域進(jìn)行精細(xì)網(wǎng)格劃分，使求解域更好地逼近各種復(fù)雜的裂縫形態(tài)，并且使用有限單元法進(jìn)行模擬計(jì)算，實(shí)現(xiàn)對(duì)線性流和徑向流酸化的模擬仿真，很好地解決了網(wǎng)格取向問(wèn)題。結(jié)果表明，將有限單元法用以酸化模擬是可行的，為多孔介質(zhì)內(nèi)的精細(xì)酸化模擬提供了新方法。

關(guān)鍵詞：裂縫；碳酸鹽巖；酸化；有限單元法；模擬計(jì)算

對(duì)于裂縫性碳酸鹽巖酸化，國(guó)內(nèi)外做了大量研究。李勇明［1］、牟建業(yè)［2］等人在進(jìn)行裂縫性碳酸鹽巖酸化模擬時(shí)，將研究重點(diǎn)放在裂縫中蚓孔的增長(zhǎng)和酸液對(duì)裂縫壁面的刻蝕。趙立強(qiáng)［3］、柳明［4］等人考慮了酸液在裂縫與地層孔隙介質(zhì)之間的相互流動(dòng)。但是大多數(shù)文獻(xiàn)中所考慮的模型都較為理想化，多是將求解域離散為規(guī)則網(wǎng)格的形式，假設(shè)裂縫與基質(zhì)均勻切割，采用有限差分法或有限體積法進(jìn)行模擬計(jì)算。然而，當(dāng)求解域中存在復(fù)雜形態(tài)的裂縫，裂縫和基巖之間幾何形態(tài)和物性參數(shù)的巨大差異以及網(wǎng)格取向性問(wèn)題，使規(guī)則網(wǎng)格形式逼近這類裂縫有較大困難。

有限元網(wǎng)格的形狀可以是任意的，常應(yīng)用于力學(xué)［5-7］和滲流問(wèn)題［8-12］，在化學(xué)動(dòng)力學(xué)問(wèn)題上也有少量應(yīng)用［13-15］，一定程度上說(shuō)明了該方法在類似領(lǐng)域的可行性。牟行洋［16］通過(guò)有限單元法對(duì)含反應(yīng)項(xiàng)的對(duì)流擴(kuò)散方程進(jìn)行了數(shù)值求解，其計(jì)算結(jié)果與解析結(jié)果吻合程度非常高，但其所求解的問(wèn)題只有數(shù)學(xué)模型而沒(méi)有結(jié)合具體實(shí)例。在前人研究的基礎(chǔ)上，使用有限元網(wǎng)格對(duì)裂縫性碳酸鹽巖區(qū)域進(jìn)行精細(xì)網(wǎng)格劃分，并使用相應(yīng)的有限單元法對(duì)其中的酸液流動(dòng)過(guò)程進(jìn)行模擬，進(jìn)行深入研究。

1　模型建立

1.1 模型假設(shè)條件

假設(shè)酸液在多孔介質(zhì)中的滲流為單相流，裂縫和基巖中的流體滲流滿足達(dá)西定律和質(zhì)量守恒定律，忽略酸巖放熱反應(yīng)造成的溫度變化對(duì)反應(yīng)速度的影響。

1.2 數(shù)學(xué)模型

多孔介質(zhì)中酸液濃度的分布受對(duì)流擴(kuò)散作用以及反應(yīng)過(guò)程的共同影響。酸溶蝕巖石后，巖石的孔隙度增大，滲透率和比面等參數(shù)也隨之變化，這些參數(shù)之間的關(guān)系往往用經(jīng)驗(yàn)公式進(jìn)行表示。Kalia［17］用雙重尺度模型來(lái)表示這一系列關(guān)系。

1.2.1 達(dá)西尺度模型 壓力場(chǎng)模型即滲流微分方程，其表達(dá)式為

式中，K為滲透率，mD；μ為酸液黏度，mPa·s；p為壓力，MPa；Ct為流體和巖石的總壓縮系數(shù)，MPa-1；為孔隙度，%；t為時(shí)間，s。

根據(jù)物質(zhì)平衡原理，某一微元內(nèi)濃度的變化等于對(duì)流擴(kuò)散產(chǎn)生的濃度變化加上反應(yīng)過(guò)程產(chǎn)生的濃度變化，其表達(dá)式為

其中

式中，Cf為酸液濃度，mol/L；u、ux、uy分別為流體流速、x方向的流速、y方向的流速，m/s；D、Dx、Dy分別為擴(kuò)散系數(shù)、x方向的擴(kuò)散系數(shù)、y方向的擴(kuò)散系數(shù)，m2/s；ks為表面反應(yīng)速率常數(shù)，m/s；kc為傳質(zhì)系數(shù)，m/s；α為酸的溶解能力，g/mol；S為比面，m-1；αos、λx和λy為與巖石結(jié)構(gòu)有關(guān)的常數(shù)；r為孔隙半徑， μm。

流速u由所得到的壓力場(chǎng)依據(jù)達(dá)西定律計(jì)算

式中，?p為壓力梯度，MPa/m。

孔隙度變化為

式中，ρ為巖石密度，kg/m3。

1.2.1 孔隙尺度模型 Kalia的孔隙尺度模型反映了滲透率、比面以及孔隙半徑等物性參數(shù)隨孔隙度的變化。滲透率響應(yīng)為

式中，K0為初始滲透率，mD；0為初始孔隙度，%；β為滲透率與孔隙度關(guān)系指數(shù)，通過(guò)實(shí)驗(yàn)獲得。

其中

式中，Sh為舍伍德數(shù)；Sh∞為漸進(jìn)舍伍德數(shù)；m為孔隙長(zhǎng)度與直徑比；Re為雷諾數(shù)；um為邊界流體流速， m/s；r0為初始孔隙半徑， μm；ν為運(yùn)動(dòng)黏度，cm2/s；Sc為施密特?cái)?shù)。

孔隙半徑與滲透率和孔隙度之間的關(guān)系式為

比面與孔隙度和孔隙半徑之間的關(guān)系式為式中，S0為初始比面，m-1。

1.3 定解條件

第一類邊界條件：流出端為恒定壓力邊界，流入端為恒定濃度邊界（t＞0）。

系統(tǒng)電源電路主要作用是為PS12034功率模塊提供 +15 V電源，為延時(shí)電路提供+24 V電源，為控制電路提供+5 V電源。選用輸入交流AC380 V/50 Hz輸出帶有+24 V、+15 V和+5 V的開(kāi)關(guān)電源。

第二類邊界條件：流入端為恒定流速邊界。

初始條件：初始時(shí)刻整個(gè)區(qū)域中的壓力為一定值，壓力濃度為0。

由于對(duì)徑向流的模擬采用的仍是以上數(shù)學(xué)模型，所以這里沒(méi)有周期邊界條件。對(duì)于封閉邊界，邊界流速為0，模型在用有限單元法離散的積分形式中的線積分項(xiàng)也為0，所以也省略了封閉邊界條件。

2　模型求解

2.1 裂縫性碳酸鹽巖物理模型

2.1.1 線性流模型的網(wǎng)格劃分 以酸液在碳酸鹽巖中的二維流動(dòng)為例，假設(shè)裂縫形狀為交叉裂縫和弧形裂縫，巖心大小為5 cm×2.5 cm，酸液從左端流入，從右端流出，類似于碳酸鹽巖巖心的酸液流動(dòng)實(shí)驗(yàn)。通過(guò)設(shè)置超大孔隙度區(qū)域（孔隙度為0.9）來(lái)模擬裂縫，設(shè)置小孔隙度區(qū)域（平均孔隙度為0.12，變化范圍為0.02~0.22）來(lái)模擬致密區(qū)域。

計(jì)算區(qū)域的網(wǎng)格劃分如圖1所示，在裂縫附近進(jìn)行了加密處理。采用三角形單元?jiǎng)澐终麄€(gè)求解域，有限元網(wǎng)格可以很好地逼近這類存在復(fù)雜裂縫的情況。雖然一般情況下，四邊形單元與三角形單元相比有更好的計(jì)算精度，但三角形單元能更好地適應(yīng)各種復(fù)雜邊界以及其他形狀較為曲折或有尖角的曲面，為了使網(wǎng)格較好地逼近各種復(fù)雜形態(tài)的裂縫，文中均使用三角形單元來(lái)進(jìn)行網(wǎng)格劃分。為清楚地看到裂縫附近的網(wǎng)格劃分情況，將圖1中紅色矩形圈出部分放大，得到局部放大圖如圖2所示。

圖1　 裂縫性碳酸鹽巖網(wǎng)格劃分

圖2　 裂縫性碳酸鹽巖網(wǎng)格劃分局部放大

2.1.2 徑向流模型的網(wǎng)格劃分 酸液由井筒向地層的流動(dòng)為徑向流，由于有限元網(wǎng)格的靈活性，這里研究酸液流動(dòng)時(shí)與Kalia、柳明等人不同，并沒(méi)有建立極坐標(biāo)系下的酸液流動(dòng)反應(yīng)模型，坐標(biāo)系仍然采用的是直角坐標(biāo)系，徑向流條件下的網(wǎng)格形式與線性流也沒(méi)有區(qū)別，井眼半徑和傷害帶半徑分別設(shè)為0.12 m和1 m，如圖3所示。

圖3　 近井多條裂縫網(wǎng)格劃分

2.2 數(shù)學(xué)模型離散

2.2.1 壓力場(chǎng)模型離散 對(duì)求解域中的任意一個(gè)單元e，首先對(duì)式（1）使用加權(quán)余量法為

其中

ωl=Nl（14）

p=Nipi+Njpj+Nmpm（15）

l=i,j,m （16）

式中，ωl為權(quán)函數(shù)；p為該單元內(nèi)的壓力近似插值函數(shù)；pi、pj、pm分別為單元節(jié)點(diǎn)的假設(shè)壓力值，MPa；Ni、Nj、Nm為基函數(shù)。

將式（13）寫(xiě)成其對(duì)應(yīng)的等效積分弱形式，并代入式（14）、（15）、（16），得到該單元內(nèi)的半離散格式為式中，Ke、Me和Fe分別為計(jì)算得到的單元?jiǎng)偠染仃?、單元溫度矩陣和單元載荷向量；Pe為未知的單元壓力向量。

按時(shí)間向后差分格式計(jì)算為

式中，?t為時(shí)間步長(zhǎng)，s；Kp、Mp和Fp分別為計(jì)算得到的總體剛度矩陣、總體溫度矩陣和總體載荷向量；Pn-1和Pn分別為上個(gè)時(shí)步和該時(shí)步的壓力場(chǎng)向量，MPa。該時(shí)步內(nèi)Pn-1為已知。

2.2.2 濃度場(chǎng)模型離散 對(duì)于濃度分布模型式（2），其中的速度u則通過(guò)所得的壓力場(chǎng)由式（5）計(jì)算，類似壓力場(chǎng)模型的離散過(guò)程，同樣可以得到

式中，Kc、Mc和Fc分別為計(jì)算濃度場(chǎng)的總體剛度矩陣、總體溫度矩陣和總體載荷向量；Cn-1和Cn分別為上個(gè)時(shí)步和該時(shí)步的濃度場(chǎng)向量。

得到求解域中的濃度分布后，則可按式（6）~（12）計(jì)算孔隙度和滲透率等物性參數(shù)，并將其代入下個(gè)時(shí)步繼續(xù)計(jì)算。

3　計(jì)算實(shí)例

根據(jù)上述模型和方法，使用Matlab編制模擬計(jì)算程序，利用表1中的數(shù)據(jù)進(jìn)行模擬計(jì)算。

表1　模型參數(shù)

3.1 線性流模擬

在相同注酸量的前提下，圖1模型分別以1×10-6m/s、1×10-5m/s和1×10-4m/s的注酸速率進(jìn)行注酸，溶解形態(tài)如圖4所示。注酸速率較低時(shí)（1×10-6m/s），端面溶蝕程度較大，酸液作用距離有限；注酸速率較高時(shí)（1×10-4m/s），雖然更多區(qū)域被溶蝕，但是無(wú)法得到足夠大的滲透率；在注酸速率適中時(shí)（1×10-5m/s），形成具有相當(dāng)導(dǎo)流能力的溶蝕通道。

為了評(píng)價(jià)酸化效果，優(yōu)選工藝參數(shù)，將巖心（或傷害帶）突破的標(biāo)志定義為入口壓力降為初始值的1%，將無(wú)因次突破體積定義為巖心突破時(shí)的酸液注入體積與巖心（或傷害帶）孔隙體積之比。如圖5所示，隨著注酸速率的增加，突破體積呈現(xiàn)先減小后增加的趨勢(shì)，即注酸速率存在一個(gè)最優(yōu)值，這里為1×10-5m/s。

3.2 徑向流模擬

圖3模型分別以1×10-5 m/s、5×10-4 m/s和5×10-3 m/s的注酸速率注酸，不同溶解形態(tài)如圖6所示，可以看出，溶解形態(tài)同樣形成了面溶蝕、蚓孔和均一溶蝕。低注酸速率（1×10-5m/s）下靠近井筒的裂縫縫寬增大，端面溶蝕程度大，但酸液作用距離有限；高注酸速率（5×10-3m/s）下形成了多分支孔道，連通多條裂縫，但不能得到足夠大的滲透率；注酸速率（5×10-4m/s）適中時(shí)形成一條高滲通道。如圖7所示，得到最優(yōu)注酸速率為5×10-4 m/s。

圖4　 線性流在不同注酸速率下的溶解形態(tài)

圖5　 線性流無(wú)因次突破體積與注酸速率曲線

圖6　 徑向流在不同注酸速率下的溶解形態(tài)

圖7　 徑向流無(wú)因次突破體積與注酸速率曲線

4　結(jié)論和建議

（1）利用有限單元法建立了與雙重尺度模型相對(duì)應(yīng)的數(shù)值模型，研究了裂縫性碳酸鹽巖儲(chǔ)層酸化特點(diǎn)。裂縫性碳酸鹽巖儲(chǔ)層酸化存在最優(yōu)注酸速率，該速率下酸液突破巖心或傷害帶所需用酸量最少。雖然原始儲(chǔ)層中存在的微裂縫幾乎不能為地下流體的滲流作貢獻(xiàn)，酸化作業(yè)可能會(huì)使這些裂縫彼此連通從而形成高滲孔道。

（2）有限元網(wǎng)格能較好地逼近各種復(fù)雜形態(tài)的裂縫，可同時(shí)實(shí)現(xiàn)對(duì)線性流和徑向流酸化的模擬，解決了網(wǎng)格取向問(wèn)題。

（3）暫時(shí)沒(méi)有考慮礦物成分分布的影響，且所研究的仍是二維模型，下一步需要開(kāi)展考慮了礦物成分的三維模型研究，做到真正意義上的仿真模擬。

參考文獻(xiàn)：

［1］ 李勇明，郭建春，趙金洲，等.裂縫性油藏酸液濾失模型研究［J］.西南石油學(xué)院學(xué)報(bào)，2004， 26（2）：50-53.

［2］ 牟建業(yè)，張士誠(chéng).酸壓裂縫導(dǎo)流能力影響因素分析［J］.油氣地質(zhì)與采收率，2011，18（2）：69-71，79.

［3］ 趙立強(qiáng)，劉平禮，劉向東，等.裂縫性碳酸鹽巖基質(zhì)酸化設(shè)計(jì)計(jì)算方法研究［J］.天然氣工業(yè)，2001，21（4）：69-71.

［4］ 柳明，張士誠(chéng)，牟建業(yè)，等.碳酸鹽巖油藏非均質(zhì)性對(duì)蚓孔擴(kuò)展的影響［J］.特種油氣藏，2012，19（5）：146-150.

［5］ 王志堅(jiān)，鄧衛(wèi)東，林忠超，等.水平井封隔器卡瓦的有限元分析及結(jié)構(gòu)改進(jìn)［J］.石油鉆采工藝，2013，35（4）：78-81.

［6］ 陳占鋒，朱衛(wèi)平，狄勤豐，等.偏磨套管破裂強(qiáng)度三維彈塑性有限元分析［J］.石油鉆采工藝，2013，35（6）：1-4.

［7］ 張輝，高德利，段明星，等.欠平衡鉆井條件下地層造斜特性研究［J］.石油鉆采工藝，2012，34（2）：1-3.

［8］ 趙敏.油藏中人工裂縫滲流規(guī)律的有限元法分析［D］.中國(guó)科學(xué)院研究生院：滲流流體力學(xué)研究所，2006.

［10］ 李亞軍，姚軍，黃朝琴，等.考慮滲透率張量的非均質(zhì)油藏有限元數(shù)值模擬方法［J］.計(jì)算物理，2010，27（5）：692-698.

［11］ 姚軍，張娜，黃朝琴，等.非均質(zhì)油藏多尺度混合有限元數(shù)值模擬方法［J］.石油學(xué)報(bào)，2012，33（3）：442-447.

［12］ 許廣明，孔祥言，盧德唐.利用非穩(wěn)定滲流方程研究射孔參數(shù)對(duì)產(chǎn)能的影響［J］.石油鉆采工藝，2000，22（4）：42-45.

［13］ KASHID M N, AGAR D W, TUREK S. CFD modelling of mass transfer with and without chemical reaction in the liquid–liquid slug flow microreactor［J］. Chemical Engineering Science, 2007, 62（18）： 5102-5109.

［14］ SEDDEEK M A. Finite-element method for the effects of chemical reaction, variable viscosity, thermophoresis and heat generation/absorption on a boundary-layer hydromagnetic flow with heat and mass transfer over a heat surface［J］. Acta Mechanica, 2005, 177（1）: 1-18.

［15］ ZHAO C, HOBBS B E, MüHLHAUS H B, et al. Finite element modelling of dissipative structures for nonequilibrium chemical reactions in fluid-saturated porous media［J］. Computer Methods in Applied Mechanics and Engineering, 2000, 184（1）: 1-14.

［16］ 牟行洋.二維對(duì)流擴(kuò)散方程的有限元求解及其反問(wèn)題研究［D］.西安：西安理工大學(xué)，2009.

［17］ KALIA N, BALAKOTAIAH V. Modeling and analysis of wormhole formation in reactive dissolution of carbonate rocks［J］. Chemical Engineering Science, 2007, 62（4）: 919-928.

（修改稿收到日期 2015-08-02）

〔編輯 李春燕〕

A new simulation method for fractured carbonate acidizing

LUO Yuan1, YANG Qianlong2, LIU Pingli1, ZHAO Liqiang1, DU Xun1, XU Bin1
（1. State Key Laboratory of Oil & Gas Reservoir Geology and Exploitation , Southwest Petroleum University, Chengdu 610500, China;
2. No.10 Oil Production Plant of Changqing Oilfield Company, CNPC, Qingyang 745100, China）

Abstract:In the previous research on simulation of fractured carbonate acidizing, it is generally assumed that the fractures cut the matrix uniformly and that the solution domain is dispersed into regular mesh pattern, which may find it very difficult to process fractures with complex shapes. In order to simulate the acid flow response process of fractured carbonate more realistically, this paper uses finite element mesh to carry out fine mesh division to its solution domain, so that the solution domain can better approach the various complex fracture shapes. In addition, finite element method is used to perform simulation, which realizes analog simulation of linear flow and radial flow acidization, favorably addressing the problem of mesh orientation. The results show that it is workable that finite element method is used for acid simulation, providing a new technique for fine acid simulation inside porous media.

Key words:fracture; carbonate; acidification; finite element method; simulation

作者簡(jiǎn)介：羅源，1987年生。西南石油大學(xué)油氣田開(kāi)發(fā)工程專業(yè)在讀碩士研究生，現(xiàn)主要從事油氣藏儲(chǔ)層改造方面的研究工作。電話：15928639637。E-mail：lyluoyuan@126.com。

基金項(xiàng)目：國(guó)家自然科學(xué)基金“復(fù)雜非均質(zhì)碳酸鹽巖儲(chǔ)層水平井酸化高效布酸基礎(chǔ)理論研究”（編號(hào)：51474182）。

doi:10.13639/j.odpt.2015.05.018

文章編號(hào)：1000 – 7393（2015）05 – 0073 – 05

文獻(xiàn)標(biāo)識(shí)碼：A

中圖分類號(hào)：TE357