考慮彌散的混溶驅替模擬

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

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

0 引言

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

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

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

1 控制方程

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

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

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

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

當圓形管束中發生混溶驅替（海水與淡水、CO2與油）時，管束中不存在完全由驅替流體所占據的部分，而變為驅替流體與被驅替流體的混合帶，此時管束的體積流量為［17］

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

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

式中：Cs為驅替流體所占據的體積分數（體積百分濃度），m3/m3。

根據節點流量守恒，流入節點的流量應等于流出節點的流量，即

根據式（3）—（5）可列出所有節點對應的線性方程組，將線性方程組轉化為矩陣的形式［式（6）］，再利用共軛梯度法對矩陣求解得到網絡模型的壓力分布

式中：A為傳導率系數矩陣；P為模型各個節點壓力組成的向量。

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

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

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

圖1 二維網格剖分示意圖Fig.1 Sketch map of 2D grid model

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

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

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

2 網絡模擬過程

基于混溶驅的基本理論及控制方程，按照以下7 個步驟進行混溶驅的網絡模擬。

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

（2）求解單相流動時網絡模型壓力場。在初始條件下，網絡模型完全飽和被驅替流體，然后結合封閉邊界條件和設定的驅替流量，求解模型的壓力場分布，此時的壓力場分布實際為被驅替流體單相流動時的壓力分布。

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

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

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

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

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

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

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

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

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

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

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

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

則式（3）也可被改寫為

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

（7）重復步驟（4），（5）和（6）實現流體的不斷驅替，當流出物中驅替流體的相對濃度C/C0達到0.999 時結束模擬。

3 模擬結果

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

3.1 不考慮和考慮彌散的混溶驅分析

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

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

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

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

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

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

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

3.2 考慮彌散的混溶驅影響因素分析

3.2.1 孔隙非均質性

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

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

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

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

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

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

3.2.2 孔隙連通性

孔隙連通性和孔隙非均質性都可以表征多孔介質微觀孔隙結構的特征，二者對混溶驅替過程的影響是類似的，但孔隙非均質性對混溶驅的影響則更大。

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

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

3.2.3 黏度比

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

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

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

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

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

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

3.2.4 驅替流量

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

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

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

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

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

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

4 結論

（1）在體中心網絡模型中通過對模型壓力場和濃度場的耦合求解，提出了考慮彌散影響的混溶驅動態網絡模擬方法。

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

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