本發(fā)明涉及一種多孔介質(zhì)中預(yù)交聯(lián)凝膠懸浮液微觀流動(dòng)的數(shù)值模擬方法，屬于油氣田開發(fā)提高原油采收率、油藏?cái)?shù)值模擬和計(jì)算流體動(dòng)力學(xué)的交叉領(lǐng)域。

背景技術(shù)：

預(yù)交聯(lián)凝膠顆粒驅(qū)是高含水油田進(jìn)一步提高原油采收率的新興深部調(diào)驅(qū)技術(shù)之一。然而，預(yù)交聯(lián)凝膠懸浮液在多孔介質(zhì)中的流動(dòng)十分復(fù)雜，不僅具備剛性顆粒懸浮液的基本滲流特征，同時(shí)還具有彈性變形等特性。目前，考慮這些特性的油藏?cái)?shù)值模擬技術(shù)多基于經(jīng)典滲濾理論或尺寸排斥理論，屬于宏觀數(shù)值模擬的范疇，受限于連續(xù)介質(zhì)假設(shè)，無法模擬預(yù)交聯(lián)凝膠顆粒在多孔介質(zhì)中“運(yùn)移、封堵、變形、再運(yùn)移”的真實(shí)過程。

為此，本發(fā)明提出一種多孔介質(zhì)中預(yù)交聯(lián)凝膠懸浮液微觀流動(dòng)的模擬方法，兼顧了顆粒變形的精細(xì)化處理、顆粒之間的接觸作用模擬以及顆粒流體間的高效耦合，為深入認(rèn)識預(yù)交聯(lián)凝膠懸浮液的復(fù)雜滲流機(jī)理提供了一種研究手段。

技術(shù)實(shí)現(xiàn)要素：

為彌補(bǔ)現(xiàn)有技術(shù)的不足，本發(fā)明將燃燒算法和轉(zhuǎn)折點(diǎn)算法結(jié)合，提出一種可以識別孔喉末端的孔隙喉道識別方法，其識別過程能準(zhǔn)確識別末端并且顯著降低識別后的喉道長度與實(shí)際長度的相對誤差。

本發(fā)明技術(shù)方案具體步驟如下：

(1)對圖像進(jìn)行二值化處理，區(qū)分巖石顆粒與孔喉空間；

(2)以巖石顆粒為前景像素對孔喉空間進(jìn)行歐式長度變換，長度變換計(jì)算公式為：

其中p為孔喉空間內(nèi)一點(diǎn)，q為p點(diǎn)周圍任意的巖石顆粒點(diǎn)。

(3)由基于長度變換的轉(zhuǎn)折點(diǎn)算法提取圖像中線，此時(shí)的中線未考慮末端；

(4)基于轉(zhuǎn)折點(diǎn)中線，進(jìn)行不考慮末端時(shí)的孔喉識別；

(5)利用燃燒算法提取中線，查找燃燒中線尾端端點(diǎn)，確定末端所在區(qū)域；

(6)借助燃燒中線尾端端點(diǎn)，利用轉(zhuǎn)折點(diǎn)算法提取末端中線；

(7)將步驟(4)中的轉(zhuǎn)折點(diǎn)中線與步驟(6)中的末端轉(zhuǎn)折點(diǎn)中線合并，對末端區(qū)域進(jìn)行孔喉識別，然后將末端識別結(jié)果與步驟(4)中未考慮末端區(qū)域的孔喉識別結(jié)果合并，得到考慮末端的孔喉識別結(jié)果。

所述步驟(4)中，孔喉識別的步驟如下：

①刪除位于孔隙內(nèi)部的轉(zhuǎn)折點(diǎn)中線路徑，將剩余中線路徑所對應(yīng)像素按順序編號，并以每個(gè)像素的長度變換結(jié)果作為此像素對應(yīng)的孔隙半徑；

②以像素編號為橫坐標(biāo)，對應(yīng)的孔隙半徑為縱坐標(biāo)，做出每條路徑的沿程曲線；

③將沿程曲線看作信號，選擇離散dmey波為小波基函數(shù)，對信號進(jìn)行二尺度分解，選擇合理閾值截取小波系數(shù)，然后重構(gòu)信號完成小波消噪(參考文獻(xiàn)：孫延奎.小波分析及其應(yīng)用[m].dynomediainc.,2005.)，閾值選擇公式為：

其中σ為噪聲標(biāo)準(zhǔn)方差，n是信號長度；

④在去噪后的沿程曲線上查找最靠近兩端的兩個(gè)局部極小值，將其對應(yīng)的像素作為孔隙和喉道的分界點(diǎn)即喉道的瓶頸，在瓶頸對應(yīng)的位置畫出孔喉識別線；

⑤孔喉識別線和巖石邊界將非基質(zhì)空間分成了若干封閉區(qū)域，將包含中線結(jié)點(diǎn)的區(qū)域作為孔隙，其余部分作為喉道。

所述步驟(6)中，提取末端中線步驟如下：

①從步驟(4)識別后的孔隙或喉道中，提取包含燃燒中線端點(diǎn)像素的喉道或孔隙區(qū)域即為末端區(qū)域；

②以燃燒中線的末端像素為原點(diǎn)，沿正交方向?qū)⑻崛〉哪┒藚^(qū)域切分為四個(gè)子圖，應(yīng)用轉(zhuǎn)折點(diǎn)算法分別提取四個(gè)子圖的中線；

③將四個(gè)子圖的中線合并在一起并清理分支，所得結(jié)果即為轉(zhuǎn)折點(diǎn)中線在末端處的分支。

本發(fā)明具有以下有益效果及優(yōu)點(diǎn)：

(1)以轉(zhuǎn)折點(diǎn)算法為主要算法提取孔隙空間的中線，使中線的整體居中效果得到改進(jìn)；

(2)利用小波消噪后的沿程曲線查找瓶頸，識別后的喉道長度與實(shí)際喉道長度相對誤差顯著降低；

(3)由燃燒算法獲取的端點(diǎn)引導(dǎo)轉(zhuǎn)折點(diǎn)算法提取末端部分的中線，可以將末端部分單獨(dú)識別出來。

附圖說明

圖1為多孔介質(zhì)中預(yù)交聯(lián)凝膠懸浮液微觀流動(dòng)的模擬方法流程圖。

圖2為多孔介質(zhì)模擬區(qū)域示意圖。

圖3為多孔介質(zhì)中數(shù)值網(wǎng)格節(jié)點(diǎn)的類型標(biāo)記示意圖。

圖4為多孔介質(zhì)中預(yù)交聯(lián)凝膠懸浮液微觀流動(dòng)的模擬結(jié)果圖。

具體實(shí)施方式

圖1給出了本發(fā)明的詳細(xì)實(shí)施流程，下面將結(jié)合其他附圖和具體實(shí)施例作進(jìn)一步的說明，以便更好的理解本發(fā)明，但不限定本發(fā)明的范圍。

(1)基于二維巖心ct切片圖像，建立多孔介質(zhì)模擬區(qū)域，如圖2所示。該模擬區(qū)域多孔介質(zhì)部分長度為260μm，為了模擬預(yù)交聯(lián)凝膠顆粒進(jìn)入和離開多孔介質(zhì)的過程，分別在多孔介質(zhì)入口端和出口端設(shè)置長度為40μm的管道，因此模擬區(qū)域總長度為340μm，寬度為100μm。設(shè)定網(wǎng)格步長為0.2μm，將該模擬區(qū)域進(jìn)行數(shù)值網(wǎng)格離散，則共有1700×500＝850000個(gè)計(jì)算節(jié)點(diǎn)。

(2)設(shè)定預(yù)交聯(lián)凝膠顆粒的初始位置在多孔介質(zhì)入口端，原始形狀為圓形，直徑為16μm。將顆粒的真實(shí)曲線邊界離散生成18個(gè)虛擬邊界標(biāo)記點(diǎn)，并由其連線組成的18邊形近似表征預(yù)交聯(lián)凝膠顆粒。

(3)標(biāo)記數(shù)值網(wǎng)格節(jié)點(diǎn)類型，若某數(shù)值網(wǎng)格節(jié)點(diǎn)及其相鄰8個(gè)節(jié)點(diǎn)均被預(yù)交聯(lián)凝膠顆粒覆蓋則為顆粒內(nèi)部節(jié)點(diǎn)(圖3中實(shí)心圓點(diǎn))，若某數(shù)值網(wǎng)格節(jié)點(diǎn)被預(yù)交聯(lián)凝膠顆粒覆蓋但存在未被覆蓋的相鄰節(jié)點(diǎn)則為顆粒邊界節(jié)點(diǎn)(圖3中實(shí)心三角形)，若某數(shù)值網(wǎng)格節(jié)點(diǎn)及其相鄰8個(gè)節(jié)點(diǎn)均未被預(yù)交聯(lián)凝膠顆粒覆蓋則為流體內(nèi)部節(jié)點(diǎn)(圖3中空心圓點(diǎn))，若某數(shù)值網(wǎng)格節(jié)點(diǎn)未被預(yù)交聯(lián)凝膠顆粒覆蓋但存在被覆蓋的相鄰節(jié)點(diǎn)則為流體邊界節(jié)點(diǎn)(圖3中空心三角形)。圖3中數(shù)值網(wǎng)格節(jié)點(diǎn)a為流體內(nèi)部節(jié)點(diǎn)，b為顆粒內(nèi)部節(jié)點(diǎn)，s為流體邊界節(jié)點(diǎn)。

(4)輸入網(wǎng)格步長δx＝0.2μm，時(shí)間步長δt＝0.0025s，松弛時(shí)間τ＝0.65，根據(jù)預(yù)交聯(lián)凝膠顆粒邊界上數(shù)值網(wǎng)格節(jié)點(diǎn)的動(dòng)量變化計(jì)算流體對各虛擬邊界標(biāo)記點(diǎn)的作用力，其計(jì)算公式為：

其中，δx表示網(wǎng)格步長，δt表示時(shí)間步長，vt表示顆粒邊界上以某數(shù)值網(wǎng)格節(jié)點(diǎn)為中心邊長為δx的控制單元的總體積，vs表示上述控制單元中被顆粒覆蓋部分的體積；τ表示松弛時(shí)間；cj為附加碰撞項(xiàng)；lj表示j方向上的離散速度；p表示各虛擬邊界標(biāo)記點(diǎn)控制范圍內(nèi)屬于顆粒邊界類型的數(shù)值網(wǎng)格節(jié)點(diǎn)個(gè)數(shù)。

(5)計(jì)算虛擬邊界標(biāo)記點(diǎn)之間的接觸重疊量δ，輸入接觸剛度k＝10⁴kg/s²，阻尼系數(shù)λ＝0.3，計(jì)算各虛擬邊界標(biāo)記點(diǎn)受到的接觸作用力，其計(jì)算公式為：

其中，k表示接觸剛度，λ表示阻尼系數(shù)，δ表示各虛擬邊界標(biāo)記點(diǎn)之間的接觸重疊量。

(6)輸入預(yù)交聯(lián)凝膠顆粒的彎曲恢復(fù)剛度γb＝10²kg/s²，拉伸恢復(fù)剛度γs＝10kg/s²，根據(jù)預(yù)交聯(lián)凝膠顆粒的當(dāng)前形狀和原始形狀，計(jì)算各虛擬邊界標(biāo)記點(diǎn)受到的彎曲恢復(fù)力和拉伸恢復(fù)力，其計(jì)算公式分別為：

fb＝γb(α-α⁰)

其中，α⁰和α分別表示預(yù)交聯(lián)凝膠顆粒變形前后相鄰三個(gè)虛擬邊界標(biāo)記點(diǎn)之間的向量夾角，和分別表示預(yù)交聯(lián)凝膠顆粒變形前后相鄰兩個(gè)虛擬邊界標(biāo)記點(diǎn)之間的距離，γb表示彎曲恢復(fù)剛度，γs表示拉伸恢復(fù)剛度。

(7)基于牛頓第二定律，計(jì)算各虛擬邊界標(biāo)記點(diǎn)的速度u和位移x，將各虛擬邊界標(biāo)記點(diǎn)相連即得到預(yù)交聯(lián)凝膠顆粒的整體形狀和位置；

x＝xold+uδt

其中，uold和u分別表示迭代前后虛擬邊界標(biāo)記點(diǎn)的速度大小，m表示預(yù)交聯(lián)凝膠顆粒的質(zhì)量，n表示虛擬邊界標(biāo)記點(diǎn)的個(gè)數(shù)，xold和x分別表示迭代前后虛擬邊界標(biāo)記點(diǎn)的位移大小。

(8)重復(fù)步驟(4)至步驟(7)，當(dāng)前后兩次模擬結(jié)果相對誤差小于0.01或者迭代次數(shù)大于100000時(shí)，結(jié)束模擬過程，輸出模擬結(jié)果，如圖4所示。

由上述實(shí)施例可以看出，本發(fā)明公開的一種多孔介質(zhì)中預(yù)交聯(lián)凝膠懸浮液微觀流動(dòng)的模擬方法能夠直觀反映預(yù)交聯(lián)凝膠顆粒在多孔介質(zhì)中的真實(shí)變形和運(yùn)移過程，并能夠定量表征預(yù)交聯(lián)凝膠懸浮液在多孔介質(zhì)中的微觀滲流機(jī)制。因此，本發(fā)明為預(yù)交聯(lián)凝膠懸浮液等含彈性顆粒懸浮液的微觀滲流機(jī)制研究提供了一種有效的數(shù)值模擬手段。