【導(dǎo)讀】注水量、開采速度等問題，因而油藏?cái)?shù)值模擬對(duì)于提高油田效益是至關(guān)緊要的。對(duì)這組方程采用有限差分方法實(shí)現(xiàn)其數(shù)值求解，以對(duì)油藏的壓力、飽和度進(jìn)行即時(shí)監(jiān)控。度采用顯式更新。經(jīng)過Matlab編寫程序，實(shí)現(xiàn)對(duì)一維油水兩相滲流方程的數(shù)值求解。

　　

【正文】 SSxQx PP nwN1nwNwv1nN1n 1Nn 1wN ????? ???? ???? ＡＡ （ 338） 于是： ???????? ????????????? ???? n1Nowwv1nN1n 1Nn 1wNn Nw1n Nw AQx PPxtSS ?? ?????? （ 339） 利用式（ 330）、式（ 332）和式（ 339）可求得 N,1,2,i ?? 網(wǎng)格的 1nwiS? 值。 以上所敘述的隱式求壓力、顯式求飽和度的計(jì)算過程即為 IMPES 方法從 n 時(shí)刻到1n? 時(shí)刻的求解過程。這樣從初始時(shí)刻開始，一步步依次求解下去，直到求得所要求的時(shí)間的壓力和飽和度為止，油水兩相的滲流方程就可得解，所建立的油藏模型就得以模擬出來。 中國(guó)石油大學(xué)（華東）本科畢業(yè)設(shè)計(jì)（論文） 19 有關(guān)單位換算 以上的計(jì)算過程中需用到達(dá)西公式，但達(dá)西公式中所用的單位為水力學(xué)單位，而實(shí)際油田中常用的單位為工程單位，二者之間的換算關(guān)系見表 31 表 31 換算關(guān)系 水力單位 工程單位 Q s/cm3 d/cm3 K ? ?Dum2 ? ?Dum2 u smPa? smPa? A 2cm 2m x? cm m P? MPa （ 1）達(dá) 西公式為xPAKQ ????，其中各變量均采用水力學(xué)單位。 （ 2）當(dāng)采用工程單位時(shí)，達(dá)西公式變?yōu)椋?43610x 10P10A10K8 64 0010Q ?????? ???。 （ 3）因此，由于水力學(xué)單位換算到工程單位，達(dá)西公式前加了一個(gè)系數(shù)，即： xPAK0 .0 8 6 4Q ??? ???? 中國(guó)石油大學(xué)（華東）本科畢業(yè)設(shè)計(jì)（論文） 20 計(jì)算程序框圖 圖 33 一維油水兩相程序框圖 程序開始 讀入 ???? 、、 xLKwo 讀入 N2,1iiKiS rww ，），（）、（ ?? 讀入 iwc PS 、 0t? 讀入 vm ax Qxt ， ? 計(jì)算不同 wS 下的 rwK 、 roK 計(jì)算 ??? 、 ow 值 計(jì)算壓力方程的系數(shù)矩陣和右邊項(xiàng) ? ? ? ? ? ? ? ?idibiaic 、 用追趕法 ??iP 顯式求 ?? ??iS,iS ow 是 打印 ? ? ? ? ? ?iSiSiPt ow 、 結(jié)束 maxtt? 否 中國(guó)石油大學(xué)（華東）本科畢業(yè)設(shè)計(jì)（論文） 21 計(jì)算實(shí)例 下面把給定的油藏區(qū)塊按一維方向均分成 7個(gè)網(wǎng)格井點(diǎn)，以第一個(gè)網(wǎng)格井為注水井，第七個(gè)網(wǎng)格井為采油井，一注一采模型進(jìn)行求解： 注水井 生產(chǎn)井 圖 34 井網(wǎng)分布圖 已知：孔隙度 ?? ，滲透率 2m1K ?? ，束縛水飽和度 ? ，原始含油飽和度 ? ，地層油粘度 smPa5o ??? ，地層水粘度 smPa1w ??? ，注水井井底壓力MPa15Piwf ? ，生產(chǎn)井井底流壓 MPa10Pwf ? ，有誰滲透率曲線如表 1 所示 .油層長(zhǎng) 70m，寬 m70 ，油層厚度 10mh? ，綜合壓縮系數(shù) 14 M P a101C ???? ， x? 取 10m ，時(shí)間步長(zhǎng)取一天，模擬時(shí)間為 1360d ，油水相對(duì)滲透率曲線見附錄。模擬曲線如下： 0 500 1000 1500 20xx00 . 20 . 40 . 60 . 81t/dfw 圖 35 最后一個(gè)網(wǎng)格點(diǎn)（生產(chǎn)井） tfw? 關(guān)系曲線 總液量中水所占有的分量即為含水率 wf ,此曲線反映了油井的含水率隨時(shí)間的變化趨勢(shì)，開始后一段時(shí)間內(nèi)含水率接近 0，因?yàn)樗?qū)前沿未到達(dá)產(chǎn)油井，采出的完全是油，水驅(qū)前沿到達(dá)產(chǎn)油井后，含水率開始上升，經(jīng)過相當(dāng)長(zhǎng)時(shí)間后含水率近似可以達(dá)到 1，此時(shí)油井開采完畢。根據(jù)此曲線可以綜合評(píng)估油井的成本與效益，以便選擇合適的開采時(shí)間。 1P 2P 3P 4P 5P 6P 7P 中國(guó)石油大學(xué)（華東）本科畢業(yè)設(shè)計(jì)（論文） 22 1 2 3 4 5 6 7111 1 . 5121 2 . 5131 3 . 5141 4 . 5iP/Mpa 圖 36 17 網(wǎng)格井點(diǎn)在 T 時(shí)刻的壓力曲線 此圖反映了此模型從第一個(gè)井點(diǎn)到第七個(gè)井點(diǎn)的壓力變化，可以看出是等差變化的。這是因?yàn)椴豢紤]油藏的毛管力 、重力以及巖石和流體的壓縮性所致，模型為剛性水驅(qū)開采，故從注水井至開采井的壓力是等差變化的。 0 200 400 600 800 1000 1200 1400 1600 18000 . 20 . 30 . 40 . 50 . 60 . 70 . 80 . 91t/dSw 圖 37 最后一個(gè)網(wǎng)格井點(diǎn)的 tSw? 曲線 此曲線反映了最后一口網(wǎng)格井的 wS 隨時(shí)間的變化趨勢(shì)，隨著水驅(qū)前沿的推進(jìn)，水侵占油藏的孔隙以排除原油，所以油藏的含水飽和度會(huì)隨時(shí)間而增大，最后達(dá)到穩(wěn)定值，因?yàn)榇藭r(shí)油 藏孔隙中的殘余油不能被排出，仍然占據(jù)著油藏孔隙，所以隨著時(shí)間的推移，含水飽和度不再變化。中國(guó)石油大學(xué)（華東）本科畢業(yè)設(shè)計(jì)（論文） 23 第 4 章 結(jié) 論 本文應(yīng)用隱式壓力 —— 顯式飽和度解法（ IMPES 方法）來對(duì)油藏油水兩相一維滲流進(jìn)行模型建立與求解，通過合并流體方程得到一個(gè)只含有壓力的方程，某一時(shí)間步的壓力求解出來后，飽和度采用顯式更新對(duì)油藏的模擬及求解。 模型的求解有定流量和定壓力兩種方法，本文建立的模型不能用定流量來進(jìn)行求解，因?yàn)槿粲枚髁窟M(jìn)行求解則此模型的壓力方程與時(shí)間無關(guān)，會(huì)有無窮多解，這是因?yàn)橹灰３肿⑺筒捎途膲毫愣?，就可?達(dá)到流量恒定，所以這樣的 壓力解是無窮的。 本文所用的 IMPES 方法具有占內(nèi)存小、計(jì)算工作量小、方法簡(jiǎn)便等優(yōu)點(diǎn)。 本文建立的模型是在不考慮巖石和流體的壓縮性以及油水的粘度、忽略重力，在這種剛性模型下對(duì)油藏的模擬準(zhǔn)確度可能有一定的影響。 IMPES 方法達(dá)西項(xiàng)的系數(shù)處理是顯式的，因此對(duì)如錐進(jìn)的問題，由于井底周圍流速高，壓差變化大，而存在較大誤差，對(duì)于強(qiáng)非線性問題的適應(yīng)性也差；飽和度的計(jì)算是顯式的，當(dāng)時(shí)間步長(zhǎng)較大時(shí)，會(huì)出現(xiàn)解的不穩(wěn)定性。 通過含水率曲線表明，相當(dāng)長(zhǎng)時(shí)間后，含水率接近 100%，此時(shí)此油帶的水 驅(qū)油開采完畢。 本文對(duì)一維油水兩相滲流進(jìn)行了數(shù)值模擬，對(duì)一維各個(gè)以劃分的網(wǎng)格在個(gè)時(shí)刻的壓力、含水飽和度、含水率等方面都有了一個(gè)全面的顯示 ， 能夠?qū)τ筒剌^為準(zhǔn)確的模擬，按此方法可以對(duì)二維、三維油水兩相以及二維、三維油氣水三相滲流進(jìn)行拓展。 中國(guó)石油大學(xué)（華東）本科畢業(yè)設(shè)計(jì)（論文） 24 致 謝 非常感謝張曉東老師在我大學(xué)的最后學(xué)習(xí)階段 —— 畢業(yè)設(shè)計(jì)階段給自己的指導(dǎo)，從最初的定題，到資料收集，到寫作、修改，到論文定稿，他給了我耐心的指導(dǎo)和無私的幫助。為了指導(dǎo)我們的畢業(yè)論文，他放棄了自己的休息時(shí)間，他的這種無私奉獻(xiàn)的敬業(yè)精神令人欽佩，在此 我向他表示我誠摯的謝意。同時(shí)，感謝所有任課老師和所有同學(xué)在這四年來給自己的指導(dǎo)和幫助，是他們教會(huì)了我專業(yè)知識(shí)，教會(huì)了我如何學(xué)習(xí)，教會(huì)了我如何做人。正是由于他們，我才能在各方面取得顯著的進(jìn)步，在此向他們表示我由衷的謝意，并祝所有的老師培養(yǎng)出越來越多的優(yōu)秀人才，桃李滿天下！ 通過這一階段的努力，我的畢業(yè)論文《基于有限差分法的油水兩相滲流方程求解》終于完成了，這意味著大學(xué)生活即將結(jié)束。在大學(xué)階段，我在學(xué)習(xí)上和思想上都受益非淺，這除了自身的努力外，與各位老師、同學(xué)和朋友的關(guān)心、支持和鼓勵(lì)是分不開的。 在本論文的寫作 過程中，我的導(dǎo)師張曉東老師傾注了大量的心血，從選題到開題報(bào)告，從寫作提綱，到一遍又一遍地指出每稿中的具體問題，嚴(yán)格把關(guān)，循循善誘，在此我表示衷心感謝。同時(shí)我還要感謝在我學(xué)習(xí)期間給我極大關(guān)心和支持的各位老師 . 寫作畢業(yè)論文是一次再系統(tǒng)學(xué)習(xí)的過程，畢業(yè)論文的完成，同樣也意味著新的學(xué)習(xí)生活的開始。我將銘記我曾是一名中國(guó)石油大學(xué)學(xué)子，在今后的工作中把惟真惟實(shí)的優(yōu)良傳統(tǒng)發(fā)揚(yáng)光大。 感謝各位專家的批評(píng)指導(dǎo) 中國(guó)石油大學(xué)（華東）本科畢業(yè)設(shè)計(jì)（論文） 25 參考文獻(xiàn) [1] 哈利德 .阿齊茲 ，安東尼 .塞特瑞 . 油藏?cái)?shù)值模擬 第一版 .北京 ：石油工業(yè)出版社， 20xx. [2] Henrik Lof， Margot Gerritsen,Marco streamline 113543,20xx. [3] 謝海兵，桓冠仁，郭尚平，等 . PEBI 網(wǎng)格二維兩相流數(shù)值模擬 .石油學(xué)報(bào)， 1999，20（ 20）： 5761. [4] 韓大匡，陳欽雷，閻存章．油藏?cái)?shù)值模擬基礎(chǔ) .北京：石油出版社， 1993． [5] 謝俊 .剩余油飽和度平面分布方法研究及應(yīng)用 .西安石油學(xué)院學(xué)報(bào)， 1998,13（ 4） :4042. [6] 張建國(guó)，杜殿發(fā)，侯健 . 油氣層滲流力學(xué)．第二版．東營(yíng)：石油大學(xué)出版社，20xx. [7] 李淑霞，谷建偉 . 油藏?cái)?shù)值模擬基礎(chǔ) 東營(yíng)：石油大學(xué)出版社 20xx. [8] [9] 鄧建中，劉之行 .計(jì)算方法 .第二版 .西安交通大學(xué)出版社 20xx. [10] 胡良劍，鄧曉君 . MATLAB 實(shí)驗(yàn) 高等教育出版社 20xx. [11] 卓金武 . MATLAB 在數(shù)學(xué)建模中的應(yīng)用 . 北京航空航天大學(xué)出版社， 20xx. 中國(guó)石油大學(xué)（華東）本科畢業(yè)設(shè)計(jì)（論文） 26 附 錄 附錄 A 表 1 油水相對(duì)滲透率曲線數(shù)據(jù)表 序號(hào) wS roK rwK 1 2 3 4 5 6 7 附錄 B 程序代碼： %******模型求解 clear all clc %************賦初值 ************* N=7。 %井?dāng)?shù) T=1800。 %天數(shù) S=100。 %面積 dx=10。 WI=5。 Pwf=10。 p0=[15。15]。 fai=。 p=zeros(N,1)。 q=zeros(N,1)。 A=zeros(N,N)。中國(guó)石油大學(xué)（華東）本科畢業(yè)設(shè)計(jì)（論文） 27 K=1。Uo=5。Uw=1。 Kro=zeros(N,T+2)。 Krw=zeros(N,T+2)。 tt=linspace(0,T,T+1)。 dert=1。 Sw=+p。 %sw0 Swtab=[ ]。 Krotab=[ 0]。 Krwtab=[0 ]。 %**************循環(huán)求解 ***************** for t=1:(T+1) Krw(:,t)=interp1(Swtab,Krwtab,Sw(:,t))。 Kro(:,t)=interp1(Swtab,Krotab,Sw(:,t))。 lamdao(:,t)=K*Kro(:,t)/Uo。 lamdaw(:,t)=K*Krw(:,t)/Uw。 lamda(:,t)=lamdao(:,t)+lamdaw(:,t)。 %**********A 矩陣賦值 *************** A(1,1)=1+WI*dx/S。 A(1,2)=1。 for i=2:(N1)