油藏数值模拟刘月田上机作业

1、实用文档中国石油大学（北京）研究生油藏数值模拟课程上机作业一油藏的储层参数分布附图所示，W1井网格位置为（4，7）, W2井位置为（7，5）。任务：独立编制计算机程序完成第五章渗流问题的模拟求解。具体要求：（一）使用显式和隐式两种方法。（二） 输出油藏投产后如下时刻的压力分布、w1井底流压和w2井产油量：初始时刻、10天、1个月、2个月、1个季度、半年、1年、2年、5年；给出压力刚好达到稳 定的时刻及其压力分布、w1井底流压和w2井产油量。（三）对照（二）中内容，比较分析显式和隐式两种方法的计算过程及结果有何不同。（四） 建议计算过程用国际单位制，输出的结果中压力用MPa产量用m/d。（五）作

2、业完成形式。要求三个电子文档：1.综合结果报告，word文档；2 显式求解方法源程序；3 隐式求解方法源程序。259222200190180185310240235228210195330290270250230205350300280259222200190180185340320290310240235228210195355335315310290270250230205325300280240210215340320290260235225355335315295275255附图 油藏的储层参数分布（1 ）厚度分布=附图中数据十50,单位：m（2）渗透率分布=附图中数据，单位：10m2

3、（3）孔隙度分布=（附图中数据X 0.02+15）%，无单位。综合结果报告一.基础数据二.显式法求解以下为运用显示求解法得到的不同时刻、不同网格压力分布(1) 初始时刻202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020W1 井井底流压 Pwf= 18.6669Mpa3W2 井产油量 Q= 93.5894m/d(2) 10 天19.394219.332319.210719.072719.01061

4、8.985019.432719.359719.200118.991118.967218.951019.508919.429819.220918.745618.915718.894819.637419.581519.440619.221119.070618.808618.38118.505818.621119.807619.778919.691719.514419.271618.828617.896718.513618.724020202019.770019.534819.237618.897218.950619.02912019.764819.572419.409419.372619.3867

5、2019.902619.817819.757519.734919.7242202020202020W1 井井底流压 Pwf= 17.4126MpaW2 井产油量 Q= 54.2204m3/d(3) 1个月19.394219.332419.210819.072719.010618.98519.432819.359719.200218.991118.967318.951119.508919.429819.22118.745718.915818.894919.637419.581619.440719.221119.070618.808718.38118.505818.621219.807719.77

6、8919.691719.514519.271618.828717.896718.513718.724120202019.770119.534919.237718.897318.950719.02912019.764819.572419.409419.372719.38672019.902619.817819.757519.734919.7242202020202020W1 井井底流压 Pwf= 17.3554MpaW2 井产油量 Q= 53.8762m3/d(4) 2个月19.323319.259719.135518.994318.92918.901319.369119.294819.1329

7、18.920918.89418.875619.45719.377219.166318.688318.855718.833819.600519.544419.401919.180419.027618.767918.351318.478418.593919.788619.759919.671419.490619.245418.803417.878318.491718.700620202019.760619.520719.221418.881118.933119.01072019.758419.56319.398719.360919.37442019.900119.813819.752619.729

8、519.7184202020202020W1 井井底流压 Pwf= 17.3552MpaW2 井产油量 Q= 53.8752m3/d(5) 1个季度19.323319.259719.135518.994318.92918.901319.369119.294819.132918.920918.89418.875619.45719.377219.166318.688318.855718.833819.600519.544419.401919.180419.027618.767918.351318.478418.593919.788619.759919.671419.490619.245418.80

9、3417.878318.491718.700620202019.760619.520719.221418.881118.933119.01072019.758419.56319.398719.360919.37442019.900119.813819.752619.729519.7184202020202020W1 井井底流压 Pwf= 17.3552Mpa3W2 井产油量 Q= 53.8752m/d(6) 半年19.323319.259719.135518.994318.92918.901319.369119.294819.132918.920918.89418.875619.45719.3

10、77219.166318.688318.855718.833819.600519.544419.401919.180419.027618.767918.351318.478418.593919.788619.759919.671419.490619.245418.803417.878318.491718.700620202019.760619.520719.221418.881118.933119.01072019.758419.56319.398719.360919.37442019.900119.813819.752619.729519.7184202020202020W1 井井底流压 P

11、wf= 17.3552MpaW2 井产油量 Q= 53.8752m3/d(7) 1 年19.323319.259719.135518.994318.92918.901319.369119.294819.132918.920918.89418.875619.45719.377219.166318.688318.855718.833819.600519.544419.401919.180419.027618.767918.351318.478418.593919.788619.759919.671419.490619.245418.803417.878318.491718.700620202019

12、.760619.520719.221418.881118.933119.01072019.758419.56319.398719.360919.37442019.900119.813819.752619.729519.7184202020202020W1 井井底流压 Pwf= 17.3552Mpa3W2 井产油量 Q= 53.8752m /d(8) 2 年19.323319.259719.135518.994318.92918.901319.369119.294819.132918.920918.89418.875619.45719.377219.166318.688318.855718.83

13、3819.600519.544419.401919.180419.027618.767918.351318.478418.593919.788619.759919.671419.490619.245418.803417.878318.491718.700620202019.760619.520719.221418.881118.933119.01072019.758419.56319.398719.360919.37442019.900119.813819.752619.729519.7184202020202020W1 井井底流压 Pwf= 17.3552Mpa3W2 井产油量 Q= 53.

14、8752m/d(9) 5 年19.323319.259719.135518.994318.92918.901319.369119.294819.132918.920918.89418.875619.45719.377219.166318.688318.855718.833819.600519.544419.401919.180419.027618.767918.351318.478418.593919.788619.759919.671419.490619.245418.803417.878318.491718.700620202019.760619.520719.221418.881118.

15、933119.01072019.758419.56319.398719.360919.37442019.900119.813819.752619.729519.7184202020202020W1 井井底流压 Pwf= 17.3552MpaW2 井产油量 Q= 53.8752m3/d(10)稳定时刻t=36d19.323319.259719.135518.994318.92918.901319.369119.294819.132918.920918.89418.875619.45719.377219.166318.688318.855718.833819.600519.544419.40191

16、9.180419.027618.767918.351318.478418.593919.788619.759919.671419.490619.245418.803417.878318.491718.700620202019.760619.520719.221418.881118.933119.01072019.758419.56319.398719.360919.37442019.900119.813819.752619.729519.7184202020202020W1 井井底流压 Pwf= 17.3552MpaW2 井产油量 Q= 53.8752m3/d二.隐式法求解以下为运用显示求解法

17、得到的不同时刻、不同网格压力分布(1) 初始时刻202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020W1 井井底流压 Pwf= 18.6669Mpa3W2 井产油量 Q= 93.5893m/d(2) 10 天19.448919.391419.278119.1519.096319.075719.483219.414519.263619.065119.0519.039519.551719.476519

18、.276218.812618.992518.978919.668919.61619.482319.273319.135518.885618.467618.606418.729419.824319.797219.715319.549119.319618.889517.959318.599718.819520202019.786619.566319.283418.954619.020319.104619.781519.601319.448219.419519.437419.909819.831219.775319.756219.74792020202020W1 井井底流压 Pwf= 17.4795

19、Mpa3W2 井产油量 Q= 55.3924m /d(3) 1个月19.341819.278819.155719.015718.951718.924819.385819.31219.15118.940218.914618.897119.470719.391219.181218.704418.87318.851819.610319.554419.412619.192119.040518.78118.362418.490118.606319.793719.765119.677119.497619.253718.812117.885518.501318.711420202019.763619.525

20、519.227418.887718.940919.01922019.760719.566719.403219.366219.38022019.901119.815419.754719.73219.7211202020202020W1 井井底流压 Pwf= 17.3713MpaW2 井产油量 Q= 54.0099m3/d(4) 2个月19.325819.262219.138218.997118.931918.904319.371419.297119.135318.923418.896618.878419.458919.379119.168318.690418.857918.83619.60181

21、9.545719.403319.181919.029318.769518.352518.479618.595119.789319.760619.672119.491519.246418.804417.87918.492618.701720202019.76119.521219.222118.881818.933919.01152019.758719.563419.399219.361419.3752019.900219.813919.752819.729819.7186202020202020W1 井井底流压 Pwf= 17.3573MpaW2 井产油量 Q= 53.8897m3/d(5) 1

22、个季度19.323619.2619.135818.994618.929318.901619.369419.29519.133118.921218.894318.87619.457219.377419.166518.688518.855918.83419.600719.544519.40219.180619.027818.768118.351418.478518.59419.788719.7619.671519.490719.245518.803517.878418.491818.700720202019.760619.520719.221518.881218.933219.01082019.7

23、58519.563119.398819.360919.37452019.900119.813819.752619.729619.7184202020202020W1 井井底流压 Pwf= 17.3555MpaW2 井产油量 Q= 53.8767m3/d(6)半年19.323319.259619.135518.994218.928918.901219.369119.294719.132818.920818.893918.875619.45719.377219.166318.688318.855618.833719.600519.544419.401919.180419.027618.767918

24、.351218.478418.593919.788619.759919.671419.490619.245318.803417.878318.491718.700620202019.760619.520719.221418.881118.933119.01072019.758419.56319.398719.360919.37442019.900119.813819.752619.729519.7184202020202020W1 井井底流压 Pwf= 17.3552Mpa3W2 井产油量 Q= 53.8750m/d(7) 1 年19.323319.259619.135518.994218.9

25、28918.901219.369119.294719.132818.920818.893918.875619.45719.377219.166318.688318.855618.833719.600519.544419.401919.180419.027618.767918.351218.478418.593919.788619.759919.671419.490619.245318.803417.878318.491718.700620202019.760619.520719.221418.881118.933119.01072019.758419.56319.398719.360919.3

26、7442019.900119.813819.752619.729519.7184202020202020W1 井井底流压 Pwf= 17.3552Mpa3W2 井产油量 Q= 53.8750m /d(8) 2 年19.323319.259619.135518.994218.928918.901219.369119.294719.132818.920818.893918.875619.45719.377219.166318.688318.855618.833719.600519.544419.401919.180419.027618.767918.351218.478418.593919.788

27、619.759919.671419.490619.245318.803417.878318.491718.700620202019.760619.520719.221418.881118.933119.01072019.758419.56319.398719.360919.37442019.900119.813819.752619.729519.7184202020202020W1 井井底流压 Pwf= 17.3552MpaW2 井产油量 Q= 53.8750m3/d(9) 5 年19.323319.259619.135518.994218.928918.901219.369119.29471

28、9.132818.920818.893918.875619.45719.377219.166318.688318.855618.833719.600519.544419.401919.180419.027618.767918.351218.478418.593919.788619.759919.671419.490619.245318.803417.878318.491718.700620202019.760619.520719.221418.881118.933119.01072019.758419.56319.398719.360919.37442019.900119.813819.752

29、619.729519.7184202020202020W1 井井底流压 Pwf= 17.3552MpaW2 井产油量 Q= 53.8750m3/d(10)稳定时刻t=121d19.323319.259619.135518.994218.928918.901219.369119.294719.132818.920818.893918.875619.45719.377219.166318.688318.855618.833719.600519.544419.401919.180419.027618.767918.351218.478418.593919.788619.759919.671419.4

30、90619.245318.803417.878318.491718.700620202019.760619.520719.221418.881118.933119.01072019.758419.56319.398719.360919.37442019.900119.813819.752619.729519.7184202020202020W1 井井底流压 Pwf= 17.3552MpaW2 井产油量 Q= 53.8750m3/d四.显隐式方法对比1、结果对比结果 时间、方法p( MPaF2 ( MPaW井底流压Pwf1 (MPaW井产量3Q2 (m/d )初始显式202018.666993.

31、5894时刻隐式202018.666993.589310天显式18.745617.896717.412654.2204隐式18.812617.959317.479555.39241个月显式18.745717.896717.355453.8762隐式18.704417.885517.371354.0099稳定显式（36 天）18.688317.878317.355253.8752时间隐式（121 天）18.688317.878317.355253.8750稳定显式18.688317.878317.355253.8752后隐式18.688317.878317.355253.8750从上表中可以看出

32、： 显式计算为有条件收敛，需要取得时间步长足够小才会收敛；而隐 式计算为无条件收敛。如在计算中输入 T=0.5d在显式计算中，经程序 运行后，会因无法收敛，产生明显错误的结果；而在隐式计算中，仍然可 以得到合理且收敛的结果。一般显式方法求解简单，稳定性差；隐式格式 求解过程复杂，稳定性好。因此一般应采用隐式计算方法，以获得较准确 的结果。 显式方法计算的压力达到稳定的时间为 36天，隐式为121天；说明在相同时间步长（并且显式收敛）时，显式解法更快的达到稳定。 两种方法都达到稳定后，1）隐式一部分地层压力分布比显式低 些（由于只保留四位有效数字，看起来不是很明显）；2）相应的稳定产量 隐式较显

33、式稍低； 以上结论可以说明，在达到同样的精度时，隐式方法计算结果更为 准确。2、计算过程对比显式差分方法：假设第n时间步的值已知，则方程中只有一个未知量，可以直接计算得到。在t=0的初值给定后，便可利用显式差分方程逐点逐层对整个 问题进行计算。求解步骤比较简单，但是稳定性较差，本问题中其稳定条 件为：取：。隐式差分方法：每个差分方程中含有五个未知量:不能独立求解，需要联立方程。其系数矩阵 A为五对角稀疏矩阵。求解过程较复杂，但稳定性较强。无条件稳定，时间间隔可以取大一些, 减少迭代次数，节省计算时间。本问题中取:显式差分求解程序P1,d,Pwf1,Q2=xscfqj(t)fun cti onP

34、i ni=20; u=5e-3;C=2e-4;Q1=30;Pwf2=15; dx=200; dy=200;dt=1.7;n=t*24/dt;re=0.208*dx;rw=0.1;K=0 259 222 200 190 180 185 0 0 0 0;259 259 222 200 190 180 185 185 0250 230 205 197.5340 320 290 3100 0;310310 240 235 228 210 195 195 0 0 0;330330 290 270180 185 0;350350 300 280 259 222 200 190 180 185 185;34

35、0240 235 228 210 195 195;355 355 335 315 310 290 270 250 230 205 205;0 00 0 325 300 280 240 210 215 215;0 0 0 0 340 320 290 260 235 225 225;0 00 0 355 335 315 295 275 255 0;H=K/50;kxd=(K.*0.02+15)/100;P=Pi ni*o nes(10,11);P1= Pi ni*on es(10,11);a=zeros(10,11);b=zeros(10,11);c=zeros(10,11);d=zeros(10

36、,11);e=zeros(10,11);for i=2:9for j=2:10a(i,j)=3600e-9*dt*2*H(i,j-1)*K(i,j-1)*K(i,j)/u/C/kxd(i,j)/dx/dx/(H(i, j-1)*K(i,j-1)+H(i,j)*K(i,j);b(i,j)=3600e-9*dt*2*H(i,j+1)*K(i,j+1)*K(i,j)/u/C/kxd(i,j)/dx/dx/(H(i, j+1)*K(i,j+1)+H(i,j)*K(i,j);c(i,j)=3600e-9*dt*2*H(i+1,j)*K(i+1,j)*K(i,j)/u/C/kxd(i,j)/dx/dx/(

37、H(i+ 1,j)*K(i+1,j)+H(i,j)*K(i,j);d(i,j)=3600e-9*dt*2*H(i-1,j)*K(i-1,j)*K(i,j)/u/C/kxd(i,j)/dx/dx/(H(i- 1,j)*K(i-1,j)+H(i,j)*K(i,j);e(i,j)=1-a(i,j)-b(i,j)-c(i,j)-d(i,j);endendtn=1;while (tn<n)dP=zeros(10,11);for i=3:6P1(i,2)=c(i,2)*P(i+1,2)+(a(i,2)+e(i,2)*P(i,2)+b(i,2)*P(i,3)+d(i,2)*P(i-1,2);dP(i,

38、2)= P1(i,2)-P(i,2);P1(i,1)=P1(i,2);endP1(2,2)=c(2,2)*P(3,2)+(a(2,2)+e(2,2)+d(2,2)*P(2,2)+b(2,2)*P(2,3);dP(2,2)= P1(2,2)-P(2,2);P1(2,1)=P1(2,2);P1(1,2)=P1(2,2);for j=3:6P1(2,j)=c(2,j)*P(3,j)+a(2,j)*P(2,j-1)+(e(2,j)+d(2,j)*P(2,j)+b(2,j)*P(2,j+1)；dP(2,j)= P1(2,j)-P(2,j);P1(1,j)=P1(2,j);endP1(2,7)=c(2,7

39、)*P(3,7)+a(2,7)*P(2,6)+(e(2,7)+b(2,7)+d(2,7)*P(2,7);dP(2,7)= P1(2,7)-P(2,7);P1(1,7)=P1(2,7);P1(2,8)=P1(2,7);for i=3:4P1(i,7)=c(i,7)*P(i+1,7)+a(i,7)*P(i,6)+(e(i,7)+b(i,7)*P(i,7)+d(i,7)*P(i-1,7)；dP(i,7)= P1(i,7)-P(i,7);P1(i,8)= P1(i,7);endfor j=8:9P1(5,j)=c(5,j)*P(6,j)+a(5,j)*P(5,j-1)+(e(5,j)+d(5,j)*P

40、(5,j)+b(5,j)*P(5,j+1)；dP(5,j)= P1(5,j)-P(5,j);P1(4,j)= P1(5,j);endP1(5,10)=c(5,10)*P(6,10)+a(5,10)*P(5,9)+(e(5,10)+b(5,10)+d(5,10)*P(5,10);dP(5,10)= P1(5,10)-P(5,10);P1(4,10)=P1(5,10);P1(5,11)=P1(5,10);for i=6:9P1(i,10)=c(i,10)*P(i+1,10)+a(i,10)*P(i,9)+(e(i,10)+b(i,10)*P(i,10)+d(i,10)*P(i-1,10);dP(i

41、,10)= P1(i,10)-P(i,10);P1(i,11)=P1(i,10);endfor i=8:9for j=6:9P1(i,j)=c(i,j)*P(i+1,j)+a(i,j)*P(i,j-1)+e(i,j)*P(i,j)+b(i,j)*P(i,j+1)+d(i,j)*P(i-1,j);dP(i,j)= P1(i,j)-P(i,j);endendfor j=5:9P1(7,j)=c(7,j)*P(8,j)+a(7,j)*P(7,j-1)+e(7,j)*P(7,j)+b(7,j)*P(7,j+1)+d (7,j)*P(6,j);dP(7,j)= P1(7,j)-P(7,j);endfor

42、 j=3:7P1(6,j)=c(6,j)*P(7,j)+a(6,j)*P(6,j-1)+e(6,j)*P(6,j)+b(6,j)*P(6,j+1)+d(6,j)*P(5,j);dP(6,j)= P1(6,j)-P(6,j);endP1(6,9)=c(6,9)*P(7,9)+a(6,9)*P(6,8)+e(6,9)*P(6,9)+b(6,9)*P(6,10)+d(6, 9)*P(5,9);dP(6,9)= P1(6,9)-P(6,9);P1(6,8)=c(6,8)*P(7,8)+a(6,8)*P(6,7)+e(6,8)*P(6,8)+b(6,8)*P(6,9)+d(6,8 )*P(5,8)-(3

43、600e-9)*2*pi*K(6,8)*dt/u/C/kxd(6,8)/dx/dy/log(re/rw)*(P(6, 8)-Pwf2);dP(6,8)= P1(6,8)-P(6,8);for j=3:7P1(5,j)=c(5,j)*P(6,j)+a(5,j)*P(5,j-1)+e(5,j)*P(5,j)+b(5,j)*P(5,j+1)+d (5,j)*P(4,j);dP(5,j)= P1(5,j)-P(5,j);endfor i=3:4for j=3:6if i=4&&j=5P1(i,j)=c(i,j)*P(i+1,j)+a(i,j)*P(i,j-1)+e(i,j)*P(i,j

44、)+b(i,j)*P(i,j+1) +d(i,j)*P(i-1,j)-dt*Q1/24/C/kxd(i,j)/H(i,j)/dx/dy;dP(i,j)= P1(i,j)-P(i,j);elseP1(i,j)=c(i,j)*P(i+1,j)+a(i,j)*P(i,j-1)+e(i,j)*P(i,j)+b(i,j)*P(i,j+1) +d(i,j)*P(i-1,j);dP(i,j)= P1(i,j)-P(i,j);endendendif max(max(abs(dP)<0.000001break ;endP=P1;tn=tn+1;endd=(t n-1)*dt/24;Pwf1=P(4,5)-

45、Q1*u*log(re/rw)/(2*pi*K(4,5)*H(4,5)/86400*1e9;Q2=2*pi*K(6,8)*H(6,8)*(P(6,8)-Pwf2 )/(u*log(re/rw)*1e-9*86400;隐式差分求解程序fun cti on wdt,p,pwf,Q=lmplict(T)dx=200;dy=200;u=5;C=2;Q 1=30;pwf2=15;K=zeros(9);H=zeros(9);fai=zeros(9);a=zeros(9);b=zeros(9);c=zeros(9);d =zeros(9);e=zeros(9);p=20* on es(9);K(1:6,4:

46、9)=355,340,350,330,310,259;335,320,300,290,240,222;315,290,2 80,270,235,200;310,310,259,250,228,190;290,240,222,230,210,180;270,235.200.205.195.185 ;K(4:6,1:3)=355,340,325;335,320,300;315,290,280;K(7:9,1:6)=295,260,240,250,228,190;275,235,210,230,210,180;255,225,215.205.195.185 ;H=K/50;fai=0.01*(0.0

47、2*K+15);HK=H.*K;re=0.208*dx;rw=0.1;dT=24;for i=1:9for j=1:9end end end endif i=9 & HK(i+1,j)=0 & HK(i,j)=0;b(i,j)=0.0864*2*HK(i+1,j)*HK(i,j)/(HK(i+1,j)+HK(i,j)/u;if i=1 & HK(i-1,j)=0 & HK(i,j)=0;a(i,j)=0.0864*2*HK(i-1,j)*HK(i,j)/(HK(i-1,j)+HK(i,j)/u;if j=9 & HK(i,j+1)=0 & HK(

48、i,j)=0;d(i,j)=0.0864*2*HK(i,j+1)*HK(i,j)/(HK(i,j+1)+HK(i,j)/u;if j=1 & HK(i,j-1)=0 & HK(i,j)=0;c(i,j)=0.0864*2*HK(i,j-1)*HK(i,j)/(HK(i,j)+HK(i,j-1)/u; g(i,j)=C*H(i,j)*dx*dy*fai(i,j)/dT*1e-4;endende=-a-b-c-d-g;J=2*3.14159*HK(7,5)/u/log(re/rw)*0.0864;nm ax=ceil(T/dT);tt=0;for n=1: nmaxnp=p；for i=1:9for j=1:9ydx(i,j)=-g(i,j)*p(i,j);endend ydx(4,7)=ydx(4,7)+Q1; ydx(7,5)=-J*pwf2+ydx(7,5);for k=1:100kp=p;for j=2:3for i=5:8p(i,j)=-(c(i,j)*p(i,j-1)+a(i,j)*p(i-1,j)+b(i,j)*p(i+1,j)+d(i,j)*p(i,j +1)-ydx(i,j)/e(i,j);endp(9,j)=-(c(9,j)*p(9,j-1)+a(9,j)*p(8,j)