一維徑向流數(shù)值模擬

1、1一維徑向單相流數(shù)學(xué)模型 對(duì)于單井問(wèn)題，通常將井底周?chē)牧鲃?dòng)看作一維徑向流， 是井底周?chē)牧髁看?、壓力變化快，而遠(yuǎn)離井底處流量小、壓力變化小，因此采 用不等距網(wǎng)格。 為模擬一維徑向單相流，首先要恰當(dāng)?shù)慕⑵鋽?shù)學(xué)模型，模型的假設(shè)如下： (1)一維徑向流動(dòng)； (2)單相流體且微可壓縮； (3)不考慮巖石的壓縮性(即巖石不可壓縮，？=常數(shù))； (4)油藏是均質(zhì)的，即 k,?為常數(shù)，流體粘度小也為一常數(shù)。 (5)不考慮重力的影響。 根據(jù)質(zhì)量守恒原理建立的柱坐標(biāo)系下單相流的數(shù)學(xué)模型為： 工鳥(niǎo)/卓)目(/需)+孰/也畢)二今“(1-1) ra尹點(diǎn)r(jOFOzzdzdt 當(dāng)只存在徑向滲流時(shí)，一維徑向單相流

2、的數(shù)學(xué)模型可簡(jiǎn)化為： 型6吟)二如以)(1-2) rc?rrcrt 考慮均質(zhì)油藏、流體微可壓縮、巖石不可壓縮，上述數(shù)學(xué)模型可簡(jiǎn)化為： 10|若用二加若用二加 c1p(1-3) rOrrdrt 假設(shè) k,?,灶勻?yàn)槌?shù)，則上述方程可簡(jiǎn)化為： 1g(孕)也空(1-4) rc，rcr，kct 方程為(1-4)即為所求的一維徑向單相流的數(shù)學(xué)模型。方程中的未知量為 p(r,t),通過(guò)求解可得沿徑向上各點(diǎn)的壓力分布及其隨時(shí)間的變化。 初始條件為:此時(shí)最典型的特點(diǎn) P(r,0)=pi(rwr0)(1-6) ar 2)定壓外邊界： p(re,t)二Pe(t0)(1-7) (2)內(nèi)邊界： 1)定產(chǎn)內(nèi)邊界： (t

3、0)(1-8) (t0)(1-9) |. (0rw21kh 2)定流壓內(nèi)邊界: P(rw,t)二Pwf 式中，r-徑向半徑，cm； rw-井底半徑，cm； re-邊界半徑,cm； p-油藏中各點(diǎn)的壓力，10-1MPa; pi-初始油藏壓力，10-1MPa; pwf-井底流壓，10-1MPa； t-時(shí)間，s； ?-孔隙度，小數(shù); k-滲透率，(jm2; C-流體的壓縮系數(shù)，1/MPa； 火流體粘度，mPa 會(huì); h-油層厚度，cm； Q-井的產(chǎn)量，cm3/s； 滲流微分方程(1-4)與初始條件、邊界條件一起，構(gòu)成了一維徑向單相流問(wèn) 題完整的數(shù)學(xué)模型。通過(guò)求解可得在各種不同的內(nèi)、外邊界條件下，地層

4、中各點(diǎn) 的壓力分布，以及井底流壓 pwf或產(chǎn)量。 2差分方程的建立 為適應(yīng)一維徑向流井底壓力變化快、遠(yuǎn)離井底附近壓力變化慢的特點(diǎn)，網(wǎng)格劃分采用不等距網(wǎng)格，即井底附近網(wǎng)格劃分密一些，遠(yuǎn)離井底要疏一些。 在此選取等比級(jí)數(shù)網(wǎng)格，即： 立二立二 a,二二 a,色二色二 a,1|1(2-1) rwri2 于是： 2._3n riarw,2araw,3abarw,%aw(2-2) 這樣實(shí)現(xiàn)了井底附近網(wǎng)格小，而遠(yuǎn)離井底處網(wǎng)格壓大的問(wèn)題。 對(duì)方程(1-4)左端項(xiàng)進(jìn)行差分，進(jìn)行一系列的變換處理，可得： Pi+iPi Pi* 1ddp、.1g1+0.50.5(4每每1-&ri)-50.5(MiAri_i)

5、 ()(G).(2-3) ror(；rrdrr 擊擊 我們將不等距的 r 坐標(biāo)轉(zhuǎn)換成了等距離的 x坐標(biāo)。兩種坐標(biāo)之間的對(duì) 應(yīng)關(guān)系如圖 1所小。則: ln乜二乜二 lna, rw 1n2.2lnaJU,rw lna=ix ln旦二灰二x1,rw lnrn.nlna,rw lnr=2Ax=x2,I,rw rni ln一二ngx二xn rw (2-4) (2-5) 于是, (2-7) r_.、一 dx1 而ln-：x為方程dx二1的特解，因此數(shù)學(xué)模型（1-4）的左端項(xiàng)可化為:rwdrir 于是數(shù)學(xué)模型（2-4）可轉(zhuǎn)換為: i62p-知Cdp 2-2-,、 C?xk0t 將式（2-8）代入上式，得:

6、1爻重 (| rdr )二二 (2-9) (2-10) r2 C?p- r2- x 22x rwie 陽(yáng)CGp I*- kdt (2-11) 通過(guò)上述過(guò)程，將不等距的徑向坐標(biāo) r 轉(zhuǎn)換成了等距離的 x坐標(biāo)，而且將數(shù) 學(xué)模型中的微分方程也進(jìn)行了坐標(biāo)轉(zhuǎn)換。 卜面用隱式差分格式對(duì)轉(zhuǎn)換為等距離 x 式(2-15)即為一維徑向流時(shí)的差分方程表達(dá)式。當(dāng) i 和？x確定以后，根據(jù)上 式用追趕法解三對(duì)角方程矩陣方程(也可直接求解)，即可確定任一半徑處的壓力分布。 3一維徑向單相流模擬事例 3.1 模擬條件與要求 已知井徑 rw=0.1m，外徑 re=250m，流體粘度.=1mPa厚度 h=5m, 滲透率 k=

7、0.05加2,孔隙度？=0.25,綜合壓縮系數(shù) C=5X10-3MPa-1,原始?jí)毫?pi=10MPa,最大模擬時(shí)間 tmax=360d,時(shí)間步？t=30d,網(wǎng)格數(shù) n=30. 外邊界定壓 p|r=re=10MPa,內(nèi)邊界定產(chǎn) Q=15m3/d。求各點(diǎn)網(wǎng)格點(diǎn)在不同時(shí)刻的壓力分布，并繪圖表示 t=90,180,270,360d時(shí)各網(wǎng)格點(diǎn)的壓力沿徑向的分布情況。 左邊的微分方程(2-11)進(jìn)行差分求解。 方程(2-11)的隱式差分方程為： pM-2pT+Ri；.22i,4x威心pn11-p re 永2.wek3t 令 22ix如CAx2 Mi：e- kAt 則式(2-12)為： n-1n11n|1

8、n pi1-(2+Mi)pipN1二-MiR (2-12) (2-13) (2-14) 令 n Zi-2Mi,di-Mipi 則: n-1；n1n|1_ pi1-/Lipipi1-di (2-15) 3.2 系數(shù)矩陣的構(gòu)建 根據(jù) 3.1 中給定的條件，可知本事例采用外邊界定壓，內(nèi)邊界定產(chǎn)的邊界條件，該類(lèi)邊界條件一般形式為： p(re,t)=Pe 逸).卜(3-1) (Sr)Lrw-2ffkhJ 下面主要構(gòu)建在上述邊界條件下，方程(2-15)對(duì)應(yīng)于 i=0 到 n 的各個(gè)網(wǎng)格所構(gòu)成的線性代數(shù)方程組。 (1)當(dāng) i=0 時(shí)，即內(nèi)邊界處，首先將內(nèi)邊界條件(r2p)|二半轉(zhuǎn)換為 x坐 r_wZikh

9、標(biāo)。轉(zhuǎn)換式如下： 0pQ“ 喘“。二而(3-2) 上式的差分方程為: Q -pwf+P1二而小 令黑-ix=do，則方程(3-3)可簡(jiǎn)化為： 2/jkh -Pwf+Pi二d0 (2)當(dāng) i=1到 n-2時(shí),按方程(2-15)列方程。 當(dāng) i=n-1 時(shí)，由式(2-15)可得： Pn-An.1Pn.1二二 dq-Pe(3-5) 當(dāng) i=n時(shí)，Pn=Pe已知，因此只需要求第 0 到 n-1 個(gè)網(wǎng)格點(diǎn)的壓力。 如上所示，列出 i=0,1,?,n-1 各網(wǎng)格節(jié)點(diǎn)的方程，所得方程組為: 當(dāng) i=0 時(shí)：Pwf+P1=do 當(dāng) i=1到 n-2時(shí)：P 售-%P1+P?；=di 當(dāng) 1 二口-1時(shí)：P.2九

10、.1Pn|1dn.1-Pe 寫(xiě)出矩陣方程的形式，得： 解此三對(duì)角矩陣方程，可求得 Pwf,P1,P2,?,Pn-1O 4計(jì)算程序框圖 一維徑向流程序框圖如圖 2所示。(3-3) (3-4) (3-6) 1 n.2 圖2一維徑向流程序框圖 5模擬結(jié)果分析 根據(jù)以上推導(dǎo)的計(jì)算公式和程序框圖，應(yīng)用 matlab進(jìn)行編程求解。主要的程序包括主程序Main、求解程序 Solve和追趕法程序 fcatch。其中，主程序 Main主要作用是輸入地層、流體參數(shù)以及初始和邊界條件，設(shè)置與模擬時(shí)間相關(guān)的參數(shù)，通過(guò)調(diào)用 Solve函數(shù)，返回一系列的結(jié)果，繪制網(wǎng)格劃分示意圖、 各網(wǎng)格點(diǎn)在不同時(shí)刻壓力分布圖和不同時(shí)刻各

11、網(wǎng)格點(diǎn)的壓力沿徑向的分布圖。 Solve函數(shù)主要作用是基于一定的邊界條件構(gòu)造系數(shù)矩陣，并調(diào)用追趕法對(duì)壓力 矩陣方程進(jìn)行求解。 為清楚顯示網(wǎng)格分布情況，根據(jù) 3中給定的條件，繪制的網(wǎng)格分布劃分示意 圖如圖 3所示。 各網(wǎng)格節(jié)點(diǎn)在不同時(shí)刻的壓力分布如圖 4所示: 圖4各網(wǎng)格點(diǎn)在不同時(shí)刻壓力分布 由圖中看出，隨網(wǎng)格編號(hào)增加（即離井越來(lái)越遠(yuǎn)），壓力下降幅度越小, 格號(hào)加 O4050 3 O3050 2 O20O 5 1 圖3網(wǎng)格劃分示意圖 小 網(wǎng)編增 I.5I,4I,3I,2 9999 10,9,876 9999 apMrp 下降的速度也越慢，在外邊界處壓力保持恒定。這是因?yàn)殡x井越遠(yuǎn)，壓力波 傳播到的

12、時(shí)間越晚，井的生產(chǎn)對(duì)該處的壓力影響也就越小 選取 t=90,180,270,360d 共 4個(gè)時(shí)間節(jié)點(diǎn)，觀察各網(wǎng)格點(diǎn)的壓力 沿徑向的分布情況。為更好的展示得到的結(jié)果，繪制 4個(gè)時(shí)間節(jié)點(diǎn)處壓力等 值線填充圖和各網(wǎng)格點(diǎn)壓力沿徑向分布圖，分別如圖 5和圖 6 所示 圖64時(shí)間節(jié)點(diǎn)各網(wǎng)格點(diǎn)壓力沿徑向分布圖 有以上兩圖中可以看出，對(duì)于圖 5,由于壓力變化較小，不同時(shí)間節(jié)點(diǎn)下的 壓力分布等值線圖相差不大，特別是當(dāng)時(shí)間為 180、270和 360d時(shí)，這一點(diǎn)在 圖 6 中體現(xiàn)的較為明顯， 三種情況下的壓力分布曲線幾乎重合。 這也說(shuō)明了生產(chǎn)之初壓力下降速t=90d t=360d 9.8 9.6 9.4 9.2 200 100 0 -100 -200 -200 0 200 9.8 9.6 9.4 9.2 圖54個(gè)時(shí)間節(jié)點(diǎn)處壓力分布等值線填充圖 050 100 150r/m 200250 300 度較快， 但由于是定壓邊界， 當(dāng)壓力波重播到邊界， 有充足的能量供給， 壓力下降速度逐漸放緩，整個(gè)油藏逐步達(dá)到穩(wěn)態(tài)。 通過(guò)以上圖像分析可知，得到的變化趨勢(shì)符合一維徑向定壓邊界油藏的開(kāi)發(fā)動(dòng)態(tài)規(guī)律，這也說(shuō)明了所編的程序是正確的，模型是合理有效的。1 上述差分格式中，由于在井底附近 ri較小，則1很大，因此易