油藏数值模拟大作业.doc

1 差分方程组线性化二维空间油水两相流体流动方程组二维空间油水两相流体流动的差分方程组油相： 水相： 其中： 令： 上式两边同乘以，则油水两相的差分方程为：油相：水相：令： 上两式可化为：差分方程线性化：方程右端项处理：油相：水相：利用IMPES方法，压力取隐式形式，饱和度和毛管压力取显示形式。令进一步可化简为：其中： 油方程/Co2水方程/Cw2，得：（1）式中：如果网格块上有井，需加上产量项，有两种情况：1.定产量，直接在左端加上Vij*qij，注入为正，生产为负;2.定井底流压，处理如下：分别将产水量和产油量加到对应的油、水方程中，即：2 求解差分方程首先，对网格进行编号，见表2.1。表2.1 网格编号510151611162145522712172226313846535759381318232732394754586049141924283340485520252934414956303542364350374451其中，18号井是注水井，11号井和46号井为生产井。将网格编号后，就可计算出公式（1）的各个系数，并将其赋值到60*60的矩阵方程中，从而可计算出各网格的压力变换量，最终求得不同生产时间的油水两相压力分布和各相饱和度分布，以及生产井的生产数据。（其中矩阵方程的系数形状见图1）图2.1 系数矩阵形状 利用C语言编制程序实现以上步骤，详细程序见附录。3 求解结果（1）分别解出100d时的各网格的油相压力及各网格的含水饱和度，见表3.1和表3.2。并根据表中数据分别绘制出地层压力分布和含水饱和度分布图，见图3.1和图3.2。表3.1 100d时油相压力分布5737.545731.54 5739.20 5793.84 5738.575737.52 5842.10 5902.66 5873.96 5848.08 5746.58 5844.99 5942.44 5850.47 6056.92 6081.26 6012.14 5953.05 5826.94 5849.02 5954.82 5957.56 5849.45 5850.15 5951.32 5936.06 6059.81 6162.36 6153.01 5957.53 5950.12 5853.40 5955.85 5956.64 5851.81 5852.36 5954.85 5958.71 6161.44 6183.26 6173.96 6060.40 6056.15 5955.60 5961.68 6162.08 6189.77 6184.64 6062.23 6058.05 5956.47 6204.12 6205.32 6104.27 6236.43 6106.21 6167.08 6257.25 6107.27 6167.54 表3.2 100d时含水饱和度分布0.4646 0.4624 0.4686 0.4672 0.4703 0.0000 0.4671 0.4739 0.4597 0.4654 0.4651 0.4663 0.4622 0.4649 0.4665 0.4683 0.4679 0.4645 0.0000 0.4671 0.4622 0.4617 0.4643 0.4622 0.4624 0.4656 0.4650 0.4656 0.4652 0.4633 0.4615 0.4634 0.4631 0.4637 0.4632 0.4629 0.4634 0.4644 0.4654 0.4651 0.4633 0.4638 0.4629 0.4638 0.4594 0.4683 0.4639 0.4644 0.4641 0.4630 0.4632 0.4655 0.4642 0.4646 0.4631 0.4639 0.4640 0.4624 0.4634 0.4630 图3.1 300d时油相压力分布 图3.2 300d时含水饱和度分布（2）分别解出300d时的各网格的油相压力及各网格的含水饱和度，见表3.3和表3.4。并根据表中数据分别绘制出地层压力分布和含水饱和度分布图，见图3.3和图3.4。表3.3 300d时油相压力分布5459.42 5450.94 5462.57 5466.24 5451.19 5436.05 5467.56 5484.22 5460.02 5480.20 5472.99 5470.78 5467.60 5480.37 5491.13 5498.93 5401.47 5488.28 5450.42 5481.75 5490.17 5494.19 5478.88 548268 5481.08 5489.16 5495.91 5400.81 5402.87 5494.87 5483.77 5488.19 5491.68 5494.31 5481.03 5482.17 5486.80 5493.59 5498.71 5402.40 5404.46 5499.31 5492.70 5491.51 5495.29 5499.87 5403.37 5405.70 5402.32 5495.60 5492.84 5404.10 5407.06 5405.86 5409.36 5409.40 5411.11 5411.10 5411.42 5412.01 表3.4 300d时含水饱和度分布0.4863 0.4493 0.5072 0.5040 0.5181 0.0000 0.5117 0.5339 0.4325 0.4858 0.4825 0.4835 0.4466 0.4801 0.4882 0.4963 0.4943 0.4888 0.0000 0.4988 0.4691 0.4645 0.4730 0.4598 0.4571 0.4800 0.4783 0.4817 0.4799 0.4709 0.4427 0.4697 0.4696 0.4731 0.4654 0.4626 0.4586 0.4731 0.4785 0.4782 0.4703 0.4709 0.4600 0.4695 0.4453 0.4897 0.4720 0.4748 0.4722 0.4605 0.4658 0.4759 0.4725 0.4729 0.4674 0.4701 0.4737 0.4642 0.4686 0.4676 图3.3 300d时油相压力分布 图3.4 300d时含水饱和度分布（3）以10天为步长，可求出生产井2和生产井3在300d内的产油量和产水量。根据这些数据绘制产量随时间的变化关系，如图5、6所示。图3.5 生产井2的产油量与产水量随时间的变化图3.6 生产井3的产油量与产水量随时间的变化附录 程序设计部分（C语言编写）#include#include#includeusing namespace std;/声明函数double Bo(double p); /原油的体积系数double Bw(double p);/水的体积系数double Oviscosity(double p);/原油的粘度double Kro(double Sw);/油的相对渗透率double Krw(double Sw);/水的相对渗透率double pc(double Sw);/毛管压力double Bo(double p) /原油的体积系数double Bo;Bo=1.0/(1.0+5.0*pow(10,-6)*(p-14.7);return Bo;double Bw(double p)/水的体积系数double Bw;Bw=1.0/(1.0+1.0*pow(10,-6)*(p-14.7);return Bw;double Oviscosity(double p)/原油的粘度double ov;int i; double ps7=5000.0,5500.0,6000.0,6500.0,7000.0,7500.0,8000.0;/压力值数组 double uo7=0.9200,0.9243,0.9372,0.9494,0.9650,0.9812,1.0019;/粘度值数组 for(i=0;i6;i+)/插值法求原油粘度 if(p=psi) ov =uoi ; else if(ppsi) continue; else ov=(p-psi)*(uoi+1-uoi)/(psi+1-psi)+uoi; break; return ov;double Kro(double Sw)/油的相对渗透率double K; int i; double swks25=0.18,0.21,0.24,0.27,0.30,0.33,0.36,0.39,0.42,0.45,0.48,0.51,0.54,0.57,0.60,0.63,0.66,0.69,0.72,0.75,0.78,0.81,0.84,0.87,0.9;/毛管力值数组 double kros25=1.00000,0.92692,0.85441,0.79288,0.71312,0.64526,0.57980,0.51709,0.45744,0.40110,0.34831,0.29924,0.25403,0.21278,0.17552,0.14228,0.11301,0.08763,0.06603,0.04803,0.3344,0.02199,0.01340,0.00733,0.0034;/油相相对渗透率值数组for(i=0;i24;i+)/插值法求油相相对渗透率 if(Sw=swksi) K=krosi ; else if(Swswksi) continue; else K=(Sw-swksi)*(krosi+1-krosi)/( swks i+1- swks i)+krosi;break; return K;double Krw(double Sw)/水的相对渗透率double K; int i; double swks25=0.18,0.21,0.24,0.27,0.30,0.33,0.36,0.39,0.42,0.45,0.48,0.51,0.54,0.57,0.60,0.63,0.66,0.69,0.72,0.75,0.78,0.81,0.84,0.87,0.9;/毛管力值数组 double krws25=0.00000,0.00000,0.00002,0.00014,0.00045,0.00111,0.00232,0.00430,0.00733,0.01175,0.01791,0.02623,0.03714,0.05116,0.06882,0.09069,0.11741,0.14963,0.18807,0.23347,0.28664,0.34842,0.41968,0.50135,0.59439;/水相相对渗透率值数组for(i=0;i24;i+)/插值法求油相相对渗透率 if(Sw=swksi) K=krwsi ; else if(Swswksi) continue; else K=(Sw-swksi)*(krwsi+1-krwsi)/( swks i+1- swks i)+krwsi;break; return K;double pc(double Sw)/毛管压力double pc; int i; double sws9=0.20,0.25,0.30,0.40,0.50,0.60,0.70,0.80,0.90;/饱和度值数组 double pcs9=8.00,4.30,3.00,1.78,1.21,0.79,0.43,0.10,0.00;/毛管力值数组 for(i=0;i8;i+)/插值法求毛管压力解 if(Sw=swsi) pc =pcsi ; else if(Swswsi) continue; else pc=(Sw-swsi)*(pcsi+1-pcsi)/(swsi+1-swsi)+pcsi; break; return pc;int main() /主函数/*/*以下是对初始变量赋初值*int i,j;double x912,y912;/数组x表示在x方向上的尺寸，数组y表示在y方向上的尺寸 double h912,D912;/数组h表示网格块的厚度，数组D表示据构造顶部深度double Kx912,Ky912;/数组Kx表示在x方向上的渗透率，数组Ky表示在y方向上的渗透率 double Porosity912;/数组Porosity表示孔隙度 double Wviscosity=0.52;/水的粘度double Cr=5.9*pow(10,-6);/综合压缩系数double p912;/原始地层压力double So912;double Sw912;double miduo=37.457;/油的密度double miduw=61.366;/水的密度double gg=32.174;/重力加速度double ro=miduo*gg*(4.42*pow(10,-5);/4.42*pow(10,-5)是单位换算系数double rw=miduw*gg*(4.42*pow(10,-5);double danwei=35.32/6.29; for(i=0;i9;i+)for(j=0;j12;j+)pij=7000.0;/原始地层压力for(i=0;i9;i+)for(j=0;j12;j+)Porosityij=0;for(i=0;i9;i+)for(j=0;j12;j+)hij=17.3;/平均厚度for(i=0;i9;i+)for(j=0;j12;j+)Dij=9307.2;/平均高度 for(i=0;i9;i+) /给x数组赋值xi0=439.0; /第一列 xi1=509.0; /第二列 xi2=491.0; /第三列 xi3=596.0; /第四列 xi4=526.0; /第五列 xi5=561.0; /第六列 xi6=912.0; /第七列 xi7=807.0; /第八列xi8=439.0; /第九列xi9=877.0; /第十列xi10=544.0; /第十一列 xi11=772.0; /第十二列 for(j=0;j12;j+) /给y数组赋值y0j=474.0;/第一行y1j=404.0; /第二行y2j=386.0; /第三行y3j=491.0; /第四行y4j=404.0; /第五行y5j=316.0; /第六行y6j=316.0; /第七行y7j=421.0; /第八行y8j=526.0; /第九行/给h数组赋值h10=8.0;h20=14.0;h30=20.0;h40=5.0;/第一列h01=10.0;h11=35.0;h21=44.0;h31=34.0;h41=12.0;/第二列h02=12.0;h12=30.0;h22=36.0;h32=35.0;h42=12.0;/第三列h03=5.0;h13=15.0;h23=30.0;h33=40.0;h43=40.0;h53=10.0;/第四列h14=6.0;h24=22.0;h34=34.0;h44=44.0;h54=19.0;/第五列h25=16.0;h35=32.0;h45=42.0;h55=27.0;h65=4.0;/第六列h26=12.0;h36=29.0;h46=32.0;h56=24.0;h66=10.0;h76=8.0;h86=4.0;/第七列h27=14.0;h37=25.0;h47=20.0;h57=10.0;h67=6.0;h77=7.0;h87=5.0;/第八列h18=4.0;h28=15.0;h38=22.0;h48=16.0;h58=6.0;h78=3.0;h88=2.0;/第九列h19=6.0;h29=11.0;h39=18.0;h49=10.0;h59=3.0;/第十列h210=6.0;h310=10.0;/第十一列h211=3.0;h311=3.0;/第十二列/给D数组赋值D10=9336.0;D20=9336.0;D30=9340.0;D40=9342.0;/第一列D01=9342.0;D11=9327.0;D21=9319.0;D31=9326.0;D41=9332.0;/第二列D02=9345.0;D12=9330.0;D22=9316.0;D32=9316.0;D42=9323.0;/第三列D03=9347.0;D13=9338.0;D23=9322.0;D33=9308.0;D43=9305.0;D53=9315.0;/第四列D14=9333.0;D24=9325.0;D34=9310.0;D44=9298.0;D54=9297.0;/第五列D25=9315.0;D35=9313.0;D45=9298.0;D55=9295.0;D65=9294.0;/第六列D26=9299.0;D36=9297.0;D46=9296.0;D56=9292.0;D66=9290.0;D76=9289.0;D86=9290.0;/第七列D27=9300.0;D37=9296.0;D47=9292.0;D57=9289.0;D67=9286.0;D77=9281.0;D87=9280.0;/第八列D18=9311.0;D28=9299.0;D38=9295.0;D48=9291.0;D58=9289.0;D78=9282.0;D88=9278.0;/第九列D19=9311.0;D29=9297.0;D39=9293.0;D49=9288.0;D59=9287.0;/第十列D210=9297.0;D310=9292.0;/第十一列D211=9305.0;D311=9295.0;/第十二列/给Kx数组赋值Kx10=267.0;Kx20=265.0;Kx30=258.0;Kx40=253.0;/第一列Kx01=275.0;Kx11=274.0;Kx21=280.0;Kx31=271.0;Kx41=259.0;/第二列Kx02=270.0;Kx12=280.0;Kx22=289.0;Kx32=295.0;Kx42=275.0;/第三列Kx03=252.0;Kx13=265.0;Kx23=278.0;Kx33=297.0;Kx43=285.0;Kx53=272.0;/第四列Kx14=253.0;Kx24=271.0;Kx34=282.0;Kx44=290.0;Kx54=276.0;/第五列Kx25=271.0;Kx35=280.0;Kx45=280.0;Kx55=273.0;Kx65=265.0;/第六列Kx26=270.0;Kx36=281.0;Kx46=289.0;Kx56=288.0;Kx66=280.0;Kx76=270.0;Kx86=260.0;/第七列Kx27=269.0;Kx37=276.0;Kx47=277.0;Kx57=281.0;Kx67=290.0;Kx77=280.0;Kx87=268.0;/第八列Kx18=259.0;Kx28=270.0;Kx38=290.0;Kx48=290.0;Kx58=274.0;Kx78=270.0;Kx88=260.0;/第九列Kx19=270.0;Kx29=279.0;Kx39=293.0;Kx49=280.0;Kx59=268.0;/第十列Kx210=283.0;Kx310=279.0;/第十一列Kx211=275.0;Kx311=270.0;/第十二列for(i=0;i9;i+)for(j=0;j12;j+)Kxij=Kxij/1000;/给Ky数组赋值Ky10=213.6;Ky20=212.0;Ky30=206.0;Ky40=202.4;/第一列Ky01=220.0;Ky11=219.2;Ky21=224.0;Ky31=216.8;Ky41=207.2;/第二列Ky02=216.0;Ky12=224.0;Ky22=231.2;Ky32=236.0;Ky42=220.0;/第三列Ky03=201.6;Ky13=212.0;Ky23=222.4;Ky33=237.6;Ky43=228.0;Ky53=217.6;/第四列Ky14=202.4;Ky24=216.8;Ky34=225.6;Ky44=232.0;Ky54=220.8;/第五列Ky25=216.8;Ky35=224.0;Ky45=224.0;Ky55=218.4;Ky65=212.0;/第六列Ky26=216.0;Ky36=224.8;Ky46=231.2;Ky56=230.4;Ky66=224.0;Ky76=216.0;Ky86=208.0;/第七列Ky27=215.2;Ky37=220.8;Ky47=221.6;Ky57=224.8;Ky67=232.0;Ky77=224.0;Ky87=214.4;/第八列Ky18=207.2;Ky28=216.0;Ky38=232.0;Ky48=232.0;Ky58=219.2;Ky78=216.0;Ky88=208.0;/第九列Ky19=216.0;Ky29=223.2;Ky39=234.4;Ky49=224.0;Ky59=214.4;/第十列Ky210=226.4;Ky310=223.2;/第十一列Ky211=220.0;Ky311=216.0;/第十二列for(i=0;i9;i+)for(j=0;j12;j+)Kyij=Kyij/1000;/单位换算/给Porosity数组赋值Porosity10=0.190;Porosity20=0.190;Porosity30=0.185;Porosity40=0.183;/第一列Porosity01=0.192;Porosity11=0.195;Porosity21=0.195;Porosity31=0.195;Porosity41=0.195;/第二列Porosity02=0.197;Porosity12=0.200;Porosity22=0.205;Porosity32=0.205;Porosity42=0.205;/第三列Porosity03=0.202;Porosity13=0.204;Porosity23=0.207;Porosity33=0.213;Porosity43=0.212;Porosity53=0.210;/第四列Porosity14=0.207;Porosity24=0.210;Porosity34=0.216;Porosity44=0.218;Porosity54=0.219;/第五列Porosity25=0.216;Porosity35=0.221;Porosity45=0.225;Porosity55=0.226;Porosity65=0.225;/第六列Porosity26=0.220;Porosity36=0.225;Porosity46=0.232;Porosity56=0.235;Porosity66=0.235;Porosity76=0.232;Porosity86=0.229;/第七列Porosity27=0.223;Porosity37=0.226;Porosity47=0.232;Porosity57=0.230;Porosity67=0.230;Porosity77=0.226;Porosity87=0.220;/第八列Porosity18=0.215;Porosity28=0.215;Porosity38=0.220;Porosity48=0.225;Porosity58=0.220;Porosity78=0.217;Porosity88=0.217;/第九列Porosity19=0.205;Porosity29=0.210;Porosity39=0.215;Porosity49=0.219;Porosity59=0.216;/第十列Porosity210=0.203;Porosity310=0.207;/第十一列Porosity211=0.200;Porosity311=0.200;/第十二列for(i=0;i9;i+)for(j=0;j12;j+)if(Porosityij=0) Soij=0;/以孔隙度为记号标记含油区域else Soij=0.50;/原始含油饱和度for(i=0;i9;i+)for(j=0;j12;j+)if(Porosityij=0) Swij=0;else Swij=0.50;/原始含水饱和度 /*以上是对初始变量赋初值*/* for(int ttt=1;ttt=300;ttt+)double t=1;/时间步长/*/*以下是求系数矩阵* double V912;double Co1912;double Co2912;double Cw1912;double Cw2912;/求V912for(i=0;i9;i+)for(j=0;j12;j+)if(Porosityij=0) Vij=0;else Vij=xij*yij*hij/5.615;/求Co1912for(i=0;i9;i+)for(j=0;j12;j+)if(Porosityij=0) Co1ij=0;else Co1ij=Vij/t*(Soij*Porosityij*(Cr/Bo(pij)+5.0*pow(10,-6);Co128=Co128+142.1*Kx28*Kro(Sw28)/(log(0.2*x28/0.25)*Oviscosity(p28)*Bo(p28);/*生产井的影响*Co112=Co112+142.1*Kx12*Kro(Sw12)/(log(0.2*x12/0.25)*Oviscosity(p12)*Bo(p12);/*生产井的影响*/求Co2912for(i=0;i9;i+)for(j=0;j12;j+)if(Porosityij=0) Co2ij=0;else Co2ij=Vij/t*Porosityij/Bo(pij);/求Cw1912for(i=0;i9;i+)for(j=0;j12;j+)if(Porosityij=0) Cw1ij=0;else Cw1ij=Vij/t*(Swij*Porosityij*(Cr/Bw(pij)+1.0*pow(10,-6);/求Cw2912for(i=0;i9;i+)for(j=0;j12;j+)if(Porosityij=0) Cw2ij=0;else Cw2ij=Vij/t*Porosityij/Bw(pij);double LDSo912;/油的流动势double LDSw912;/水的流动势double Kro1912;/用上游权法判断的x方向上油的相对渗透率double Kro2912;/用上游权法判断的y方向上油的相对渗透率double Krw1912;/用上游权法判断的x方向上水的相对渗透率 double Krw2912;/用上游权法判断的y方向上水的相对渗透率double Kx1912;/用调和平均求得的x方向上绝对渗透率double Ky1912;/用调和平均求得的y方向上绝对渗透率 /油的流动势 for(i=0;i9;i+)for(j=0;j12;j+)if(Porosityij=0) LDSoij=0;else LDSoij=pij-ro*Dij;/水的流动势 for(i=0;i9;i+)for(j=0;j12;j+)if(Porosityij=0) LDSwij=0;else LDSwij=pij-rw*Dij;/用调和平均求得的x方向上绝对渗透率for(i=0;i9;i+)for(j=0;j12;j+)if(Porosityij=0)Kx1ij=0;else if(j=11) Kx1ij=Kxij; else Kx1ij=2*Kxij*Kxij+1/(Kxij+Kxij+1);/用调和平均求得的y方向上绝对渗透率 for(i=0;i9;i+)for(j=0;j12;j+)if(Porosi