考慮彌散的混溶驅(qū)替模擬

戚 濤，胡 勇，李 騫，趙梓寒，張 春，李 滔

（1.中國石油西南油氣田分公司勘探開發(fā)研究院，成都 610041；2.中國石油西南油氣田分公司，成都 610051）

0 引言

混溶驅(qū)替指2 種或2 種以上可混溶流體的驅(qū)替，既可以是單相流體（海水與淡水、地表污水與地下水、含示蹤劑的流體和純流體），也可以是兩相流體（CO2與原油等），其廣泛應(yīng)用于水文學(xué)、化學(xué)、醫(yī)學(xué)和石油工程等領(lǐng)域［1-4］，而混相驅(qū)替主要應(yīng)用于石油工業(yè)領(lǐng)域，且更側(cè)重于2 種不同相態(tài)的可混溶流體的驅(qū)替［3-7］。由于多孔介質(zhì)孔隙結(jié)構(gòu)特征（孔隙非均質(zhì)性、孔隙連通性、迂曲度等）和流體流動特征（被驅(qū)替流體與驅(qū)替流體的黏度比、分子擴散、彌散、驅(qū)替流量等）的差異，混溶驅(qū)替過程常常表現(xiàn)出黏性指進現(xiàn)象。因此，掌握影響混溶驅(qū)替的主控因素對于弄清溶質(zhì)運移規(guī)律和擴大驅(qū)替流體的體積波及效率具有重要意義。

國內(nèi)外學(xué)者［8-16］利用連續(xù)模型［8］、隨機模型［9］、孔隙網(wǎng)絡(luò)模擬［10-11］和Lattice-Boltzmann 模擬［12-14］等方式開展了大量混溶驅(qū)替研究：Mohiuddin 等［8］利用Navier-Stokes 方程、連續(xù)性方程和對流擴散方程研究了重力對混溶驅(qū)替過程的影響，該方法很難考慮多孔介質(zhì)的孔隙結(jié)構(gòu)特征，僅能進行定性分析，不能作定量評價；Tartakovsky 等［15］用粒子表示流體，運用光滑粒子流體力學(xué)方法（SPH）求解Lan‐gevin 方程，實現(xiàn)了混溶過程中對流和擴散的有效分離。Ebrahimi 等［9］運用小波變換法粗化非均勻多孔介質(zhì)來實現(xiàn)混溶驅(qū)替的升尺度模擬，發(fā)現(xiàn)滲透率各向異性和黏度比均是影響油田尺度混溶驅(qū)替效率的關(guān)鍵因素；Xia［12］重點研究了孔隙結(jié)構(gòu)、Le 數(shù)和Re 數(shù)等對混溶驅(qū)替過程的影響；Stevenson 等［10］和Sideiqui 等［16］利用孔隙網(wǎng)絡(luò)模型開展了流體的隨機性或確定性運動模擬，研究了非均質(zhì)性、流度比等因素對混溶驅(qū)的影響，但均未考慮彌散的作用。

為加深對混溶驅(qū)替微觀機理的認識，通過對網(wǎng)絡(luò)模型壓力場和濃度場的耦合求解，開展考慮彌散的混溶驅(qū)替動態(tài)網(wǎng)絡(luò)模擬，定量分析孔隙非均質(zhì)性、孔隙連通性、黏度比和驅(qū)替流量等對巖心尺度下混溶驅(qū)替的影響，以期弄清混溶驅(qū)替過程的主控因素。

1 控制方程

當圓形管束中發(fā)生單相流體流動時，體積流量服從Hagen-Poiseuille 方程［17］，即

式中：qij為管束流量，m3/s；pi，pj分別為由管束所連接的節(jié)點i，j的壓力，MPa；rij為管束半徑，m；μ為流體黏度，mPa·s；l為管束長度，m。

當圓形管束中發(fā)生不混溶驅(qū)替（油和水）時，不同流體間存在明顯的分界面，若驅(qū)替流體占據(jù)的管束長度為Xij，則被驅(qū)替流體占據(jù)的長度為l-Xij，此時管束的體積流量為

式中：μ1，μ2分別為驅(qū)替流體、被驅(qū)替流體的流體黏度，mPa·s；Xij為驅(qū)替流體或混合帶所占據(jù)的管束長度，m；。

當圓形管束中發(fā)生混溶驅(qū)替（海水與淡水、CO2與油）時，管束中不存在完全由驅(qū)替流體所占據(jù)的部分，而變?yōu)轵?qū)替流體與被驅(qū)替流體的混合帶，此時管束的體積流量為［17］

式中：μe為混合帶的流體黏度，mPa·s。

混合帶的黏度μe根據(jù)流體混合的四次方根理論［18］進行計算

式中：Cs為驅(qū)替流體所占據(jù)的體積分數(shù)（體積百分濃度），m3/m3。

根據(jù)節(jié)點流量守恒，流入節(jié)點的流量應(yīng)等于流出節(jié)點的流量，即

根據(jù)式（3）—（5）可列出所有節(jié)點對應(yīng)的線性方程組，將線性方程組轉(zhuǎn)化為矩陣的形式［式（6）］，再利用共軛梯度法對矩陣求解得到網(wǎng)絡(luò)模型的壓力分布

式中：A為傳導(dǎo)率系數(shù)矩陣；P為模型各個節(jié)點壓力組成的向量。

在確定網(wǎng)絡(luò)模型的初始壓力分布后，需要對網(wǎng)絡(luò)模型的濃度場分布進行求解。驅(qū)替流體在單個孔隙的流動滿足Taylor-Aris 彌散［19-20］且服從一維對流彌散方程［21-22］（ADE）

式中：C為體積百分濃度，m3/m3；t為時間，s；D為Taylor-Aris 彌散系數(shù)，m2/s；vm為管束的平均流速，m/s；Dm為分子擴散系數(shù)，m2/s。

采用顯式差分格式求解非線性方程，計算管束中的濃度分布。兼顧模型求解的準確性和計算的可行性，將單根管束剖分為4 個網(wǎng)格（以二維網(wǎng)絡(luò)模型為例，圖1），則空間步長△x等于管束長度的1/4。

圖1 二維網(wǎng)格剖分示意圖Fig.1 Sketch map of 2D grid model

圖1 中，若驅(qū)替流體從節(jié)點1（i，j）流入節(jié)點5（i+1，j），則管束中節(jié)點1 的濃度即可得知，而管束中的節(jié)點2，3，4 和5 對應(yīng)的濃度C2，C3，C4和C5分別為

式中：Ci（i=1，2，3，4，5）為節(jié)點i的體積百分濃度，m3/m3；?x為空間步長，m；?t為時間步長，s；上標n為第n次迭代。

顯式差分若要收斂和穩(wěn)定必須對時間步長和空間步長加以限制。訥正［23］的研究發(fā)現(xiàn)一維對流彌散方程顯式差分必須滿足式（13）和式（14），本次研究所取參數(shù)均滿足該條件。

2 網(wǎng)絡(luò)模擬過程

基于混溶驅(qū)的基本理論及控制方程，按照以下7 個步驟進行混溶驅(qū)的網(wǎng)絡(luò)模擬。

（1）構(gòu)建多孔介質(zhì)網(wǎng)絡(luò)模型。按照多孔介質(zhì)網(wǎng)絡(luò)模型的建模方法［22］構(gòu)建模型大小為40×40×5的體中心（BCC）網(wǎng)絡(luò)模型，管束半徑服從對數(shù)均勻分布，不同的配位數(shù)通常通過隨機斷開管束的方式獲得。不同參數(shù)的網(wǎng)絡(luò)模型表征不同孔隙結(jié)構(gòu)的多孔介質(zhì)。

（2）求解單相流動時網(wǎng)絡(luò)模型壓力場。在初始條件下，網(wǎng)絡(luò)模型完全飽和被驅(qū)替流體，然后結(jié)合封閉邊界條件和設(shè)定的驅(qū)替流量，求解模型的壓力場分布，此時的壓力場分布實際為被驅(qū)替流體單相流動時的壓力分布。

（3）確定計算時間步長?T。根據(jù)被驅(qū)替流體單相流動的壓力場分布計算流體通過單根管束的平均流體速度vij，Taylor-Aris 彌散系數(shù)D和平均時間tij。

式中：vij為通過單根管束ij的平均流體速度，m/s；Sij為單根管束ij橫截面積，m2；tij為通過單根管束ij的平均時間，s。

找出驅(qū)替流體通過單根管束的最短時間

式中：Tmin為驅(qū)替流體通過單根管束的最短時間，s。

以最短時間的1/10 為每次壓力更新的時間步長（或計算時間步長）?T

式中：?T為計算時間步長，s。

計算時間步長?T受模型非均質(zhì)性、連通性以及驅(qū)替流量等因素的共同影響，變化范圍大致為10-4～10-2s。由于計算時間步長?T較短，可認為在?T內(nèi)的驅(qū)替過程屬于穩(wěn)定驅(qū)替。差分時間步長?t取為計算時間步長的1/100，即10-6～10-4s，以保證顯式差分格式收斂且穩(wěn)定。

（4）更新單管中的流體分布。不考慮彌散影響時，驅(qū)替流體在?T內(nèi)向前推進的距離為

驅(qū)替流體在單管中占據(jù)的總長度為所有?Xij的疊加，即Xij=∑?Xij。若，則須對計算時間步長?T進行調(diào)整，?T取為驅(qū)替前緣剛好到達管束節(jié)點的最短時間，即

考慮彌散影響時，在所有與驅(qū)替流體接觸的管束中，用式（9）～（12）的顯式差分格式求解單根管束中各節(jié)點的濃度分布，結(jié)合流入節(jié)點的管束流量，利用加權(quán)平均的方法計算網(wǎng)絡(luò)模型節(jié)點的濃度。

（5）更新各管束傳導(dǎo)率，為下一次壓力求解作準備。不考慮彌散時，根據(jù)更新后驅(qū)替流體和被驅(qū)替流體所占管束長度更新管束傳導(dǎo)率?？紤]彌散時，單根管束為驅(qū)替流體和被驅(qū)替流體的混合帶。由圖1 可以發(fā)現(xiàn)，單根管束被剖分為4 個網(wǎng)格，共使用了5 個節(jié)點描述管束的濃度分布，每個網(wǎng)格的體積百分濃度可認為是網(wǎng)格兩端點的平均值，因此，每一個網(wǎng)格的有效黏度為

式中：μei為網(wǎng)格i的有效黏度，mPa·s；μi，μi+1分別為網(wǎng)格i，i+1 的流體黏度，mPa·s。

則式（3）也可被改寫為

（6）根據(jù)修正后的有效黏度及傳導(dǎo)率，利用共軛梯度法求解式（6）以得到定流量q所對應(yīng)的壓力場分布，隨后更新各孔隙的平均流體速度和Taylor-Aris 彌散系數(shù)Dij。

（7）重復(fù)步驟（4），（5）和（6）實現(xiàn)流體的不斷驅(qū)替，當流出物中驅(qū)替流體的相對濃度C/C0達到0.999 時結(jié)束模擬。

3 模擬結(jié)果

選取特定的網(wǎng)絡(luò)模型參數(shù)（水力半徑R=4 μm、管束長度l=200 μm、變異系數(shù)C V=0.8、配位數(shù)z=4.8）和流體運動參數(shù)（驅(qū)替流量q=10-5mL/s、被驅(qū)替流體與驅(qū)替流體的黏度比或流度比M=10、分子擴散系數(shù)Dm=10-10m2/s）作為混溶驅(qū)的標準條件，以此進行單因素分析。

3.1 不考慮和考慮彌散的混溶驅(qū)分析

彌散控制著驅(qū)替流體和被驅(qū)替流體的混合程度，不考慮彌散的混溶驅(qū)只存在對流作用，該過程與加入表面活性劑的油水非混溶驅(qū)類似，驅(qū)替過程沒有考慮毛管壓力的影響。

由圖2 可看出：當被驅(qū)替流體與驅(qū)替流體的黏度比為1 時，壓差和壓力場分布幾乎不發(fā)生變化，當黏度比大于1 時，壓差隨著注入孔隙體積（PV數(shù)）的增加而減??；在驅(qū)替剛開始時，考慮彌散的壓差下降速度更慢，這有利于保持地層能量，擴大體積波及效率，提高采收率；當注入PV 數(shù)達到0.8時，在不考慮彌散的情形中，驅(qū)替流體沿著貫穿的通道直接流出模型，被驅(qū)替流體很難再被驅(qū)替出來，壓差趨于穩(wěn)定，在考慮彌散的情形中，仍有部分被驅(qū)替流體被驅(qū)替出來，壓差進一步減小；隨著注入PV 數(shù)的進一步增加，2 種情形的壓差都趨于一個定值，但是考慮彌散時對應(yīng)的壓差更小，主要是因為穩(wěn)定情形類似于驅(qū)替流體的單相流動，考慮彌散的模型中驅(qū)替流體占據(jù)的體積相對更大，且驅(qū)替流體的黏度小于被驅(qū)替流體，導(dǎo)致最終壓差偏小。

圖2 考慮彌散和不考慮彌散時的壓差與注入PV 數(shù)的關(guān)系Fig.2 Relationship between pressure difference and injected PV number with and without dispersion

不考慮彌散時，驅(qū)替流體主要進入主流通道中［24］，表現(xiàn)出黏性指進現(xiàn)象［25］［圖3（a）］?？紤]彌散時，驅(qū)替流體以主流通道為基準向四周擴散，進入了更多的孔隙，其中包括部分流量較小的孔隙，從而導(dǎo)致沿主流通道流動的流體減少［圖3（b）］，雖然主流通道中的流體也會在彌散作用下發(fā)生向前的擴展，但是在該流量（v=10-5mL/s）的驅(qū)替下，其作用十分有限，并不能消除主流通道中流體減少帶來的影響，因此驅(qū)替流體較為均勻的向前推進，壓制了黏性指進現(xiàn)象。對比圖3（a）和圖3（b）可以發(fā)現(xiàn)，2 種情形均是從同一根管束突破，說明多孔介質(zhì)的無序性決定了流體流動的主流通道，而彌散作用只是將通道進行延展。

圖3 不考慮彌散（a）和考慮彌散（b）前緣突破時的情形對比Fig.3 Comparison of breakthrough with dispersion（a）and without dispersion（b）

圖4 為流出物中驅(qū)替流體相對濃度C/C0隨注入PV 數(shù)的變化關(guān)系，該曲線類似于油水非混溶驅(qū)的含水率曲線。由圖4 可以發(fā)現(xiàn)：不考慮彌散時，驅(qū)替流體在0.26 PV 發(fā)生突破，對應(yīng)的采出程度為26%，突破后，C/C0迅速上升，達到0.7 左右，然后再緩慢上升；考慮彌散時，驅(qū)替流體在0.36 PV 發(fā)生突破，對應(yīng)的采出程度為36%，比不考慮彌散對應(yīng)的突破時間要晚，說明彌散作用有效壓制了黏性指進現(xiàn)象，且C/C0呈緩慢上升的形態(tài)，不存在驟變的情況，當注入PV 數(shù)達到2 時，C/C0基本達到1。

圖4 C/C0 隨注入PV 數(shù)的變化Fig.4 Relationship betweenC/C0 and injected PV number

3.2 考慮彌散的混溶驅(qū)影響因素分析

3.2.1 孔隙非均質(zhì)性

圖5 為標準條件下其他參數(shù)不變時，不同變異系數(shù)下驅(qū)替前緣突破時的二維圖形。由圖5 可以看出，C為0.05 時，驅(qū)替流體均勻推進，且過渡帶很短；隨著C V增加，孔隙半徑差異性增大，黏性指進現(xiàn)象越明顯，形成的指進通道越多；隨著C V進一步增加，模型中會出現(xiàn)難以被驅(qū)替的管束，驅(qū)替流體在此類管束附近發(fā)生繞流，呈現(xiàn)中心管束未被驅(qū)替而四周管束被驅(qū)替的情形，形成殘余區(qū)域。

圖5 不同變異系數(shù)下前緣突破時的分布Fig.5 Fluid distribution of breakthrough under different coefficient of variation

由圖6 可以看出：當C V=0.05 時，前緣突破時間最晚，突破后C/C0上升的速度最快，當注入PV 數(shù)達到1.2 時，C/C0幾乎不再發(fā)生變化，且等于1；隨著C V的增加，前緣突破時間提前，而C/C0達到1 的時間延后，這大大增加了前緣突破后的生產(chǎn)時間；注入PV 數(shù)為1 所對應(yīng)的C/C0隨C V的增加稍有增加。

圖6 不同變異系數(shù)下C/C0隨注入PV 數(shù)的變化Fig.6 Relationship betweenC/C0 and injected PV number under different coefficient of variation

由圖7 可知：采出程度與C V滿足線性關(guān)系，C V為0.05 時對應(yīng)的采出程度最高，為85%；隨著C V增加，采出程度急劇下降，當C V達到1.05 時，采出程度僅為12%。

圖7 不同變異系數(shù)下前緣突破時的采出程度Fig.7 Recovery degree of breakthrough under different coefficient of variation

3.2.2 孔隙連通性

孔隙連通性和孔隙非均質(zhì)性都可以表征多孔介質(zhì)微觀孔隙結(jié)構(gòu)的特征，二者對混溶驅(qū)替過程的影響是類似的，但孔隙非均質(zhì)性對混溶驅(qū)的影響則更大。

總體來看，前緣突破時的采出程度與配位數(shù)（z）呈線性關(guān)系（圖8）。當模型完全連通（z=8.0）時，前緣突破時間最晚，C/C0到達1 的時間最早，對應(yīng)的采出程度最高，約為61%；隨著z的減小，形成的指進通道增多，殘余區(qū)域面積增大，前緣突破時間逐漸縮短，而C/C0到達1 的時間逐漸延長。

圖8 不同配位數(shù)下前緣突破時的采出程度Fig.8 Recovery degree of breakthrough under different coordination number

3.2.3 黏度比

被驅(qū)替流體和驅(qū)替流體的黏度比嚴重影響著混溶驅(qū)替過程。當黏度比（M）為2 時，整個驅(qū)替過程基本呈活塞式驅(qū)替；隨著M的增大，主流通道效應(yīng)顯著，整個驅(qū)替過程表現(xiàn)出黏性指進現(xiàn)象，且M越大，黏性指進現(xiàn)象越明顯；當M達到100 時，驅(qū)替流體幾乎只沿一條通道流動，體積波及效率極低（圖9）。

圖9 不同黏度比下前緣突破時的分布Fig.9 Fluid distribution of breakthrough under different viscosity ratio

圖10 不同黏度比下C/C0 隨注入PV 數(shù)的變化Fig.10 Relationship betweenC/C0 and injected PV number under different viscosity ratio

隨著M的增加，前緣突破時間逐漸縮短，而C/C0到達1 的時間逐漸延長（圖10）。當M增大到一定程度時（M≥50），前緣突破時間幾乎不再發(fā)生變化，但前緣突破后C/C0的上升速度和幅度受M影響，對于M為50 的情形，C/C0隨著注入PV 數(shù)的增加緩慢增加，對于M為100 的情形，C/C0快速上升到0.6 左右，然后隨著注入PV 數(shù)的增加緩慢增加。

由圖11 可看出：M越接近1，前緣突破時的采出程度越高；隨著M的增加，采出程度急劇下降，其變化趨勢基本服從乘冪規(guī)律；當M大于等于50 時，前緣突破時的采出程度開始趨于平穩(wěn)。由此可以發(fā)現(xiàn)，當M較大時，驅(qū)替流體會快速突破，使體積波及效率偏低，因此在CO2驅(qū)油過程中常加入泡沫以增加流體在高滲通道中的滲流阻力或加入高分子量聚合物以增加CO2黏度，以此防止CO2快速突破。

圖11 不同黏度比下前緣突破時的采出程度Fig.11 Recovery degree of breakthrough under different viscosity ratio

3.2.4 驅(qū)替流量

驅(qū)替流量直接決定對流作用和彌散作用的強弱。當驅(qū)替流量較小時，對流作用相對較弱，彌散作用相對較強，驅(qū)替流體在彌散作用下可大范圍向四周擴散，從而使驅(qū)替流體較為均勻地向前驅(qū)替，呈現(xiàn)活塞式驅(qū)替。隨著驅(qū)替流量的增加，對流作用大大增強，而彌散作用稍有增強，導(dǎo)致驅(qū)替流體主要在主流通道流動，向四周擴散的范圍有限，表現(xiàn)出黏性指進現(xiàn)象，且存在許多殘余區(qū)域。當驅(qū)替流量增加到一定值（模擬條件下10-4mL/s）時，驅(qū)替過程幾乎不再發(fā)生變化（圖12）。

圖12 不同驅(qū)替流量下前緣突破時的分布Fig.12 Fluid distribution of breakthrough under different displacement velocity

在模擬條件下，驅(qū)替流量小于等于10-6mL/s時，驅(qū)替過程差異很小，前緣突破時間和C/C0到達1 的時間幾乎一致，說明驅(qū)替流體受對流作用的影響十分有限，且流體在彌散作用下已經(jīng)完成了所有的橫向擴散，體積波及效率已經(jīng)達到極致。隨著驅(qū)替流量的增大，前緣突破時間提前，C/C0到達1 的時間延后。當驅(qū)替流量增大到一定程度時，對流作用將占主體地位，驅(qū)替流量對前緣突破時間的影響減弱，轉(zhuǎn)為影響最終采收率（圖13）。

圖13 不同驅(qū)替流量下C/C0隨注入PV 數(shù)的變化Fig.13 Relationship betweenC/C0and injected PV number under different displacement velocity

當注入流量小于等于10-6mL/s 時，前緣突破時的采出程度基本保持在60%不變；當注入流量大于10-6mL/s 時，采出程度隨著注入流量的增加而減小，其變化趨勢基本服從乘冪規(guī)律（圖14）。

圖14 不同驅(qū)替流量下前緣突破時的采出程度Fig.14 Recovery degree of breakthrough under different displacement velocity

4 結(jié)論

（1）在體中心網(wǎng)絡(luò)模型中通過對模型壓力場和濃度場的耦合求解，提出了考慮彌散影響的混溶驅(qū)動態(tài)網(wǎng)絡(luò)模擬方法。

（2）多孔介質(zhì)的無序性決定流體流動的主流通道，而彌散作用將主流通道進行延展。與不考慮彌散的混溶驅(qū)相比，考慮彌散的混溶驅(qū)的模型兩端壓差減小速度更慢，最終穩(wěn)定時的壓差更小，體積波及效率更大，驅(qū)替流體突破時間更晚，前緣突破時的采出程度更高，流出物中驅(qū)替流體的相對濃度C/C0上升更為緩慢。

（3）對于考慮彌散的混溶驅(qū)替過程，孔隙非均質(zhì)性、孔隙連通性、黏度比和驅(qū)替流量等對混溶驅(qū)均具有較大影響，其中黏度比對混溶驅(qū)的影響相對較強，孔隙連通性對混溶驅(qū)的影響相對較弱?？紫斗蔷|(zhì)性越強（或孔隙連通性越低或黏度比越大或驅(qū)替流量越大），前緣突破時間越早，黏性指進現(xiàn)象越明顯，體積波及效率越低。前緣突破時的采出程度與孔隙非均質(zhì)性（或孔隙連通性）呈線性關(guān)系，與黏度比（或驅(qū)替流量）呈乘冪關(guān)系。