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

1、1 一維徑向單相流數(shù)學(xué)模型對于單井問題，通常將井底周圍的流動看作一維徑向流，此時最典型的特點是井底周圍的流量大、壓力變化快，而遠(yuǎn)離井底處流量小、壓力變化小，因此采用不等距網(wǎng)格。為模擬一維徑向單相流，首先要恰當(dāng)?shù)慕⑵鋽?shù)學(xué)模型，模型的假設(shè)如下：（1） 一維徑向流動；（2） 單相流體且微可壓縮；（3） 不考慮巖石的壓縮性（即巖石不可壓縮，=常數(shù)）；（4） 油藏是均質(zhì)的，即k，為常數(shù)，流體粘度也為一常數(shù)。（5） 不考慮重力的影響。根據(jù)質(zhì)量守恒原理建立的柱坐標(biāo)系下單相流的數(shù)學(xué)模型為： (1-1)當(dāng)只存在徑向滲流時，一維徑向單相流的數(shù)學(xué)模型可簡化為：(1-2)考慮均質(zhì)油藏、流體微可壓縮、巖石不可壓縮，上

2、述數(shù)學(xué)模型可簡化為：(1-3)假設(shè)k，均為常數(shù)，則上述方程可簡化為：(1-4)方程為(1-4)即為所求的一維徑向單相流的數(shù)學(xué)模型。方程中的未知量為p(r,t)，通過求解可得沿徑向上各點的壓力分布及其隨時間的變化。初始條件為：P(r,0)=pi (rwrre)(1-5)邊界條件包括外邊界和內(nèi)邊界。相應(yīng)的外邊界條件如下：(1) 外邊界：1) 封閉外邊界：(1-6)2) 定壓外邊界： (1-7)(2) 內(nèi)邊界：1) 定產(chǎn)內(nèi)邊界：(1-8)2) 定流壓內(nèi)邊界：(1-9)式中，r-徑向半徑，cm；rw-井底半徑，cm；re-邊界半徑，cm；p-油藏中各點的壓力，10-1MPa；pi-初始油藏壓力，10-

3、1MPa；pwf-井底流壓，10-1MPa；t-時間，s；-孔隙度，小數(shù)；k-滲透率，m2；C-流體的壓縮系數(shù)，1/MPa；-流體粘度，mPas；h-油層厚度，cm；Q-井的產(chǎn)量，cm3/s；滲流微分方程(1-4)與初始條件、邊界條件一起，構(gòu)成了一維徑向單相流問題完整的數(shù)學(xué)模型。通過求解可得在各種不同的內(nèi)、外邊界條件下，地層中各點的壓力分布，以及井底流壓pwf或產(chǎn)量。2 差分方程的建立為適應(yīng)一維徑向流井底壓力變化快、遠(yuǎn)離井底附近壓力變化慢的特點，網(wǎng)格劃分采用不等距網(wǎng)格，即井底附近網(wǎng)格劃分密一些，遠(yuǎn)離井底要疏一些。在此選取等比級數(shù)網(wǎng)格，即： (2-1)于是： (2-2)這樣實現(xiàn)了井底附近網(wǎng)格小，

4、而遠(yuǎn)離井底處網(wǎng)格壓大的問題。對方程(1-4)左端項進(jìn)行差分，進(jìn)行一系列的變換處理，可得： (2-3)上述差分格式中，由于在井底附近ri較小，則很大，因此易造成計算的不穩(wěn)定，故應(yīng)將空間坐標(biāo)做適當(dāng)?shù)淖儞Q，即將一維的徑向坐標(biāo)轉(zhuǎn)換為直角坐標(biāo)。為把一維徑向坐標(biāo)r轉(zhuǎn)換為直角坐標(biāo)x，需要找到r與x的對應(yīng)關(guān)系。由式(2-2)可得： (2-4)令 則：(2-5)于是，我們將不等距的r坐標(biāo)轉(zhuǎn)換成了等距離的x坐標(biāo)。兩種坐標(biāo)之間的對應(yīng)關(guān)系如圖1所示。圖1 不等距r坐標(biāo)與等距x坐標(biāo)之間的轉(zhuǎn)換已知rw，re和網(wǎng)格數(shù)n時，可以求出轉(zhuǎn)換后的網(wǎng)格大小x。由可得： (2-6)由式(2-5)可看出，r與x之間的對應(yīng)關(guān)系為：(2-7

5、)于是： (2-8)而為方程的特解，因此數(shù)學(xué)模型(1-4)的左端項可化為： (2-9)于是數(shù)學(xué)模型(2-4)可轉(zhuǎn)換為：(2-10)將式(2-8)代入上式，得：(2-11)通過上述過程，將不等距的徑向坐標(biāo)r轉(zhuǎn)換成了等距離的x坐標(biāo)，而且將數(shù)學(xué)模型中的微分方程也進(jìn)行了坐標(biāo)轉(zhuǎn)換。下面用隱式差分格式對轉(zhuǎn)換為等距離x左邊的微分方程(2-11)進(jìn)行差分求解。方程(2-11)的隱式差分方程為： (2-12)令(2-13)則式(2-12)為： (2-14)令則：(2-15)式(2-15)即為一維徑向流時的差分方程表達(dá)式。當(dāng)i和x確定以后，根據(jù)上式用追趕法解三對角方程矩陣方程（也可直接求解），即可確定任一半徑處的

6、壓力分布。3 一維徑向單相流模擬事例3.1 模擬條件與要求已知井徑rw=0.1m，外徑re=250m，流體粘度=1mPas，厚度h=5m，滲透率k=0.05m2，孔隙度=0.25，綜合壓縮系數(shù)C=5×10-3MPa-1，原始壓力pi=10MPa，最大模擬時間tmax=360d，時間步t=30d，網(wǎng)格數(shù)n=30.外邊界定壓p|r=re=10MPa，內(nèi)邊界定產(chǎn)Q=15m3/d。求各點網(wǎng)格點在不同時刻的壓力分布，并繪圖表示t=90，180，270，360d時各網(wǎng)格點的壓力沿徑向的分布情況。3.2 系數(shù)矩陣的構(gòu)建根據(jù)3.1中給定的條件，可知本事例采用外邊界定壓，內(nèi)邊界定產(chǎn)的邊界條件，該類邊界

7、條件一般形式為：(3-1)下面主要構(gòu)建在上述邊界條件下，方程(2-15)對應(yīng)于i=0到n的各個網(wǎng)格所構(gòu)成的線性代數(shù)方程組。(1) 當(dāng)i=0時，即內(nèi)邊界處，首先將內(nèi)邊界條件轉(zhuǎn)換為x坐標(biāo)。轉(zhuǎn)換式如下： (3-2)上式的差分方程為： (3-3)令,則方程(3-3)可簡化為：(3-4)(2) 當(dāng)i=1到n-2時，按方程(2-15)列方程。(3) 當(dāng)i=n-1時，由式（2-15）可得： (3-5)(4) 當(dāng)i=n時，pn=pe已知，因此只需要求第0到n-1個網(wǎng)格點的壓力。如上所示，列出i=0，1，n-1各網(wǎng)格節(jié)點的方程，所得方程組為：當(dāng)i=0時：當(dāng)i=1到n-2時：當(dāng)i=n-1時：寫出矩陣方程的形式，得

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

9、進(jìn)行求解。為清楚顯示網(wǎng)格分布情況，根據(jù)3中給定的條件，繪制的網(wǎng)格分布劃分示意圖如圖3所示。圖3 網(wǎng)格劃分示意圖各網(wǎng)格節(jié)點在不同時刻的壓力分布如圖4所示：網(wǎng)格編號增加圖4 各網(wǎng)格點在不同時刻壓力分布由圖中看出，隨網(wǎng)格編號增加（即離井越來越遠(yuǎn)），壓力下降幅度越小，下降的速度也越慢，在外邊界處壓力保持恒定。這是因為離井越遠(yuǎn)，壓力波傳播到的時間越晚，井的生產(chǎn)對該處的壓力影響也就越小。選取t=90，180，270，360d共4個時間節(jié)點，觀察各網(wǎng)格點的壓力沿徑向的分布情況。為更好的展示得到的結(jié)果，繪制4個時間節(jié)點處壓力等值線填充圖和各網(wǎng)格點壓力沿徑向分布圖，分別如圖5和圖6所示。圖5 4個時間節(jié)點處壓力分布等值線填充圖圖6 4時間節(jié)點各網(wǎng)格點壓力沿徑向分布圖有以上兩圖中可以看出，對于圖5，由于壓力變化較小，不同時間節(jié)點下的壓力分布等值線圖相差不大，特別是當(dāng)時間為180、270和360d時，這一點在圖6中體現(xiàn)的較為明顯，