一維徑向流數值模擬

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

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

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

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

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

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

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

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

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