資源描述:
《油藏數(shù)值模擬實驗報告.doc》由會員上傳分享,免費在線閱讀,更多相關(guān)內(nèi)容在教育資源-天天文庫。
1、數(shù)值模擬上機實驗1:三對角系數(shù)矩陣的解法隱式差分格式出發(fā)點在(i,n+1),取關(guān)于t的一階向后差商和關(guān)于x的二階差商。取對于一維滲流問題的隱式差分方程組的系數(shù)矩陣為三對角矩陣,追趕法(THOMAS)就是用來求解三對角矩陣方程組的一種比較簡單、應(yīng)用也極為廣泛的解法。它的基本思路是將三對角矩陣A分解成兩個特定形式的三對角矩陣的乘積。追趕法程序如下:DimmAsInteger,nAsInteger,iAsIntegerDimP(1To10)AsSingle,x(1To10)AsSingle,y(1To10)AsSingle,a(1To10)AsSingle,b(
2、1To10)AsSingle,c(1To10)AsSingle,d(1To10)AsSingleDiml(1To10)AsSingle,u(1To10)AsSinglen=InputBox("請輸入方程個數(shù)")Fori=1Tona(i)=InputBox("a("&i&")=?")Print"a("&i&")=";a(i);NextiPrintFori=1Ton-1b(i)=InputBox("b("&i&")=?")Print"b("&i&")=";b(i);NextiPrintFori=2Tonc(i)=InputBox("c("&i&")=?")Pr
3、int"c("&i&")=";c(i);NextiPrintFori=1Tond(i)=InputBox("d("&i&")=?")Print"d("&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)Print"y(1)=";y(1);Fori=2Tony(i)=(d(i)-c(i)*y(i-1))/l(i)Print"y("&i&")=";y(i);NextiPrintx(n)=y(n)Fori=n-1T
4、o1Step-1x(i)=y(i)-u(i)*x(i+1)NextiFori=1TonPrint"x("&i&")=";x(i);Nexti運行結(jié)果:實驗2:單相流數(shù)值模擬方法已知一維均質(zhì)等厚單相滲流,已知邊界條件定壓、初始壓力分布,求任意時刻的壓力分布。滲流微分方程為:網(wǎng)格系統(tǒng)為均勻的塊中心網(wǎng)格,求t=任意時刻的壓力分布(隱式差分格式)。一維隱式求解:DimP(0To100,0To100)AsSingle,rAsSingle,nAsInteger,kAsIntegerDimiAsInteger,jAsInteger,tAsDouble,zAsSingle
5、,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=1
6、Tonc(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(
7、i)*x(i+1)NextiFori=1TonP(i,k)=x(i)NextiNextkFori=1TonPrintP(i,tm);NextiPrint運行結(jié)果:二維程序如下:DimP(0To4,0To3)AsSingle,r(1To7,1To7)AsSingle,kAsIntegerDimiAsInteger,jAsInteger,qAsInteger,tAsSingle,zAsSingle,tmAsIntegerDimx(1To6)AsSingle,y(1To6)AsSingle,rmAsIntegertm=InputBox("輸入所求時間")r(1,
8、1)=-6:r(1,2)=1:r(1,3)=1r(2,1)=1:r