油藏数值模拟实验报告

数值模拟上机实验1：三对角系数矩阵的解法隐式差分格式出发点在（i,n＋1)，取关于t的一阶向后差商和关于x的二阶差商。取对于一维渗流问题的隐式差分方程组的系数矩阵为三对角矩阵，追赶法（THOMAS）就是用来求解三对角矩阵方程组的一种比较简单、应用也极为广泛的解法。它的基本思路是将三对角矩阵A分解成两个特定形式的三对角矩阵的乘积。iniiniiniidpbpapc1111122ppxttppxpppninininini121111122txninininipppp1111121追赶法程序如下:DimmAsInteger,nAsInteger,iAsIntegerDimP(1To10)AsSingle,x(1To10)AsSingle,y(1To10)AsSingle,a(1To10)AsSingle,b(1To10)AsSingle,c(1To10)AsSingle,d(1To10)AsSingleDiml(1To10)AsSingle,u(1To10)AsSinglen=InputBox(请输入方程个数)Fori=1Tona(i)=InputBox(a(&i&)=?)Printa(&i&)=;a(i);NextiPrintFori=1Ton-1b(i)=InputBox(b(&i&)=?)Printb(&i&)=;b(i);1d2dxNdidacb1p2pxNpipNextiPrintFori=2Tonc(i)=InputBox(c(&i&)=?)Printc(&i&)=;c(i);NextiPrintFori=1Tond(i)=InputBox(d(&i&)=?)Printd(&i&)=;d(i);NextiPrintl(1)=a(1)Fori=2Tonu(i-1)=b(i-1)/l(i-1)l(i)=a(i)-c(i)*u(i-1)Nextiy(1)=d(1)/l(1)Printy(1)=;y(1);Fori=2Tony(i)=(d(i)-c(i)*y(i-1))/l(i)Printy(&i&)=;y(i);NextiPrintx(n)=y(n)Fori=n-1To1Step-1x(i)=y(i)-u(i)*x(i+1)NextiFori=1TonPrintx(&i&)=;x(i);Nexti运行结果:实验2:单相流数值模拟方法已知一维均质等厚单相渗流，已知边界条件定压、初始压力分布，求任意时刻的压力分布。渗流微分方程为：网格系统为均匀的块中心网格，求t=任意时刻的压力分布（隐式差分格式）。一维隐式求解：DimP(0To100,0To100)AsSingle,rAsSingle,nAsInteger,24122),(5),0()1(10)0,()0(22LtxtLptpLxxpLxtpxp，，取ninininipppp1111121iniiniiniidpbpapc111111d2dxNdidacb1p2pxNpipkAsIntegerDimiAsInteger,jAsInteger,tAsDouble,zAsSingle,xmAsSingle,tmAsIntegerDimx(1To100)AsSingle,y(1To100)AsSingle,a(1To100)AsSingle,b(1To100)AsSingle,c(1To100)AsSingle,d(1To100)AsSingleDiml(1To100)AsSingle,u(1To100)AsSinglez=1:xm=2:n=4r=z/xmv=-(1+2*r)tm=InputBox(输入所求时间)Fori=1To4P(i,0)=10NextiFori=0TotmP(0,i)=5:P(5,i)=2NextiFori=1Tonc(i)=r:a(i)=v:b(i)=rNextiFork=1Totmd(1)=-P(1,k-1)-b(1)*P(0,k)d(2)=-P(2,k-1)d(3)=-P(3,k-1)d(4)=-P(4,k-1)-b(1)*P(5,k)l(1)=a(1)Fori=2Tonu(i-1)=b(i-1)/l(i-1)l(i)=a(i)-c(i)*u(i-1)Nextiy(1)=d(1)/l(1)Fori=2Tony(i)=(d(i)-c(i)*y(i-1))/l(i)NextiPrintx(n)=y(n)Fori=n-1To1Step-1x(i)=y(i)-u(i)*x(i+1)NextiFori=1TonP(i,k)=x(i)NextiNextkFori=1TonPrintP(i,tm);NextiPrint运行结果：二维程序如下：DimP(0To4,0To3)AsSingle,r(1To7,1To7)AsSingle,kAsIntegerDimiAsInteger,jAsInteger,qAsInteger,tAsSingle,zAsSingle,tmAsIntegerDimx(1To6)AsSingle,y(1To6)AsSingle,rmAsIntegertm=InputBox(输入所求时间)r(1,1)=-6:r(1,2)=1:r(1,3)=1r(2,1)=1:r(2,2)=-6:r(2,4)=1r(3,1)=1:r(3,3)=-6:r(3,4)=1:r(3,5)=1r(4,2)=1:r(4,3)=1:r(4,4)=-6:r(4,6)=1r(5,3)=1:r(5,5)=-6:r(5,6)=1r(6,4)=1:r(6,5)=1:r(6,6)=-6Forq=1TotmP(1,1)=7:P(2,1)=6:P(3,1)=5P(1,2)=6:P(2,2)=5:P(3,2)=4y(1)=-2*P(1,1)-10-8y(2)=-2*P(1,2)-10-5y(3)=-2*P(2,1)-7y(4)=-2*P(2,2)-4y(5)=-2*P(3,1)-4-6y(6)=-2*P(3,2)-3-2Fori=1To6r(i,7)=y(i)NextiFori=1To6Forj=1To7Printr(i,j);NextjPrintNextiFork=1To5Fori=k+1To6t=r(i,k)/r(k,k)Forj=k+1To7r(i,j)=r(i,j)-r(k,j)*tNextjNextiNextkFori=1To6Forj=1To7Printr(i,j);NextjPrintNextiFori=6To1Step-1Forj=i+1To6z=x(j)*r(i,j)+zNextjx(i)=(r(i,7)-z)/r(i,i)NextiFori=1To3P(i,1)=x(i)P(i,2)=x(i+3)NextiNextqForj=1To2Fori=1To3PrintP(i,j);NextiPrintNextj程序结果：