《油藏數(shù)值模擬》第三章

*1油藏數(shù)值模擬NumericalReservoirSimulation主講：石油與天然氣工程學(xué)院*2第三章微分方程的離散化*3開采過程非線性偏微分方程單/多相流公式非線性代數(shù)方程離散化線性代數(shù)方程線性化①建立數(shù)學(xué)模型②建立數(shù)值模型A、通過離散化將偏微分方程組轉(zhuǎn)換為有限差分方程組。B、將非線性系數(shù)線形化，得到線形代數(shù)方程組。Ａ、通過質(zhì)量/能量守恒方程、狀態(tài)方程、運動方程、輔助方程建立基本方程組。B、根據(jù)所研究的具體問題建立相應(yīng)的初始和邊界條件。*4第三章微分方程的離散化第一節(jié)偏導(dǎo)數(shù)的差商逼近第二節(jié)網(wǎng)格系統(tǒng)第三節(jié)有限差分方程的建立--空間離散第四節(jié)差分方程的建立--時間離散第五節(jié)差分方程的穩(wěn)定性分析第六節(jié)線性代數(shù)方程組的求解方法*5典型流動方程—黑油模型*6黑油方程的典型表示忽略重力時間：一階偏微分空間：二階偏微分*7為什么要進行離散化處理？簡單的微分方程，在簡單的條件下，可以用解析法求解。實際油藏流動中，飽和度、壓力等隨時間和空間的變化不能用簡單函數(shù)來描述：油藏分布的非均質(zhì)性；相滲曲線、毛管壓力隨飽和度的關(guān)系一般是非線性的；流體的體積、組分的物理性質(zhì)和壓力、溫度的關(guān)系是非線性的。*8考察函數(shù)u(x),其自變量的一階導(dǎo)數(shù)*9考察函數(shù)u(x),其自變量的一階導(dǎo)數(shù)可定義為下面的各種極限前差商中心差商后差商一、一階偏導(dǎo)數(shù)的差商逼近第一節(jié)偏導(dǎo)數(shù)的差商逼近*101.一階向前差商一階向前差商對函數(shù)u(x+x)在x處進行Talor級數(shù)展開得:忽略一階小量0(x)得下式*112.一階向后差商一階向后差商對函數(shù)u(x-x)在x處進行Talor級數(shù)展開得:忽略一階小量0(x)得下式*123.一階中心差商一階中心差商忽略二階小量0(x2)*13二、二階偏導(dǎo)數(shù)的差商逼近忽略二階小量0(x2)標(biāo)準(zhǔn)二階差商逼近不等距網(wǎng)格上二階差商示意圖形如:的差商逼近x-x’xx+x’’(x’+x’’)/2

x’

x’’*15其中:二、二階偏導(dǎo)數(shù)的差商逼近上式即為

的差商表達式,可以證明其誤差也是取

,且，則上式變?yōu)?*17第二節(jié)網(wǎng)格系統(tǒng)油藏油藏有效網(wǎng)格無效網(wǎng)格(i,j)(i,j)ab。

。

。*18常規(guī)網(wǎng)格系統(tǒng)非常規(guī)網(wǎng)格系統(tǒng)塊中心網(wǎng)格點中心網(wǎng)格局部網(wǎng)格加密混合網(wǎng)格加密局部正交網(wǎng)格（PEBI）第二節(jié)網(wǎng)格系統(tǒng)多邊形網(wǎng)格*19一、網(wǎng)格系統(tǒng)--塊中心網(wǎng)格i,ji+1,ji-1,ji,j-1i,j+1xyxixi+1xi+1yj+1yjyj+1定流量或封閉

特點：模擬計算使用的體積與實際體積一致，且感覺也更直觀，易于被大家接受。*20二、網(wǎng)格系統(tǒng)--點中心網(wǎng)格i,ji+1,ji-1,ji,j-1i,j+1xyxixi-1xi+1yj+1yjyj-1定壓邊界

特點：邊界正好與計算點重合，因此在第一類邊界條件下有利于提高模擬的精度。*21三、非常規(guī)網(wǎng)格系統(tǒng)--局部網(wǎng)格加密*22四、非常規(guī)網(wǎng)格系統(tǒng)--混合網(wǎng)格加密*23一、一維網(wǎng)格的近似第三節(jié)有限差分方程的建立--空間離散1234x*24展開得代數(shù)方程組:一、一維網(wǎng)格的近似上式中:*25一、一維網(wǎng)格的近似*26一、一維網(wǎng)格的近似矩陣AP=f*27aibici一維自然排列時三對角矩陣結(jié)構(gòu)示意圖一、一維網(wǎng)格的近似*28第三節(jié)有限差分方程的建立--空間離散二、二維網(wǎng)格的近似五點差分格式i-1,ji,ji+1,ji,j-1i,j+1*31網(wǎng)格排序按行的自然排序按列的自然排序二、二維網(wǎng)格的近似按行的自然排序代數(shù)方程組二、二維網(wǎng)格的近似*33矩陣AP=fwi二、二維網(wǎng)格的近似aijbijcijdijeij二維自然排列時五對角矩陣結(jié)構(gòu)示意圖*35三、三維網(wǎng)格的近似類似于二維展開得代數(shù)方程組:第三節(jié)有限差分方程的建立--空間離散aijbijcijdijeij三維自然排列時七對角矩陣結(jié)構(gòu)示意圖gijfij*37四、空間離散的物理意義簡記第三節(jié)有限差分方程的建立--空間離散*38四、空間離散的物理意義累積項的變化=流入–

流出*39第四節(jié)差分方程的建立--時間離散時間離散把研究的實踐域分成若干小的時間段，在每個時間段內(nèi)，對問題求解，時間段之間有機連接。步長大小取決于所要解決的實際問題。*40第四節(jié)差分方程的建立--時間離散一、顯示差分格式將三維單相微可壓縮流體的數(shù)學(xué)模型簡化為一維拋物方程：利用P(x，t)關(guān)于t的一階向前差商和關(guān)于x的二階差商，在點（i，n）的差分方程。令

則上式變?yōu)?又因：*42從方程可以看出：如果已知第n（本步時間）的值Pin，就可以求的第n+1時刻（下步時間）的值Pin+1。因此，如初始條件，即n=0時刻各網(wǎng)格的P值已給定，就可以依次求的以后各時間的P值。在其中只有一個未知數(shù)Pin+1，由一個方程就可以求出。簡單，精度較差，時間步長受到嚴(yán)格限制，基本不用。第四節(jié)差分方程的建立--時間離散*43第四節(jié)差分方程的建立--時間離散二、隱式差分格式將三維單相微可壓縮流體的數(shù)學(xué)模型簡化為一維拋物方程：利用P(x，t)關(guān)于t的一階向后差商和關(guān)于x的二階差商，在點（i，n+1）的差分方程。*44第四節(jié)差分方程的建立--時間離散從方程可以看出：如果已知第n（本步時間）的值Pin，為了求得第n+1時刻（下步時間）的值Pin+1，必須解一個線性代數(shù)方程組。即：要想求出Pin+1值，需用到第n時刻的P值，也要用到n+1時刻的P值。這種差分格式就是隱式差分格式。在點（i，n+1），用到點（i-1，n+1)、（i+1，n+1)和（i，n)三個點。只有一個方程，卻有2-3個未知數(shù)。穩(wěn)定，精度好，廣泛使用。*45三、CN差分格式Crank_Nicolson差分格式（簡稱CN差分格式）是綜合顯式和隱式格式而構(gòu)建，將空間二階差商取為n時刻與n+1時刻的算術(shù)平均值，則有：第四節(jié)差分方程的建立--時間離散求解精度高，工作量與隱式格式差不多，在油藏數(shù)值模擬中經(jīng)常被采用。*46四、其它差分格式第四節(jié)差分方程的建立--時間離散*47五、二維不穩(wěn)定滲流方程的時間離散第四節(jié)差分方程的建立--時間離散顯式差分格式利用P(x，y，t)關(guān)于t的一階向前差商和關(guān)于x，y的二階差商，在點（i，j，n）的差分方程。*48五、二維不穩(wěn)定滲流方程的時間離散第四節(jié)差分方程的建立--時間離散該線性代數(shù)方程組在節(jié)點（i，j）列方程式，用到（i，j）,(i+1,j),(i-1,j),(i,j+1),(i,j-1)五個點：顯式：只有一個方程，1個未知數(shù)。簡單，精度較差，時間步長受到嚴(yán)格限制，基本不用。*49五、二維不穩(wěn)定滲流方程的時間離散第四節(jié)差分方程的建立--時間離散隱式差分格式利用P(x，y，t)關(guān)于t的一階向后差商和關(guān)于x，y的二階差商，在點（i，j，n+1）的差分方程。*50五、二維不穩(wěn)定滲流方程的時間離散第四節(jié)差分方程的建立--時間離散該線性代數(shù)方程組在節(jié)點（i，j）列方程式，用到（i，j）,(i+1,j),(i-1,j),(i,j+1),(i,j-1)五個點，但時刻不同。隱式：一個方程，5個未知數(shù)，穩(wěn)定且廣泛使用。*51第四節(jié)差分方程的建立--時間離散嘗試推導(dǎo)三維空間不穩(wěn)定滲流的顯式和隱式差分格式通式？*52第四節(jié)差分方程的建立--時間離散六、邊界條件的處理（一）內(nèi)邊界條件的處理定產(chǎn)條件：即井以一定產(chǎn)量q生產(chǎn)。如在網(wǎng)格（i，j）上有一口井，產(chǎn)量為q，則可在滲流方程左邊加上產(chǎn)量項，生產(chǎn)井q為負(fù)，注水井q為正。定壓條件：即井以一定流動壓力Pwf生產(chǎn)，這時的q未知，可由給定的井底流動壓力Pwf和井點所在網(wǎng)格節(jié)點的壓力Pij計算。把網(wǎng)格內(nèi)井近似堪稱穩(wěn)態(tài)流動，符合平面徑向流：？*53第四節(jié)差分方程的建立--時間離散關(guān)于供給半徑re，不同的學(xué)者，有不同的公式：1）對于各向同性地層，其等值供給半徑：2）對于各向異性地層，其等值供給半徑：*54第四節(jié)差分方程的建立--時間離散六、邊界條件的處理（二）外邊界條件的處理封閉邊界：常取塊中心網(wǎng)格并在邊界網(wǎng)格外虛擬一排網(wǎng)格，并令其相鄰兩個網(wǎng)格壓力相等。定壓邊界：常取點中心網(wǎng)格，由于邊界點的壓力一定，因此只需求內(nèi)部節(jié)點壓力。*55第五節(jié)差分方程的穩(wěn)定性分析前面計算Pi是分時刻進行的，如果再n時刻引入誤差，必然會影響到n+1步時刻以及以后各時刻的數(shù)值。一個重要的問題是：當(dāng)本步帶來誤差時，利用它去計算后面的結(jié)果，誤差將會怎樣傳播？是逐漸消失，還是逐步積累，越來越大？如果誤差越來越大，差分方法計算的解與微分方程的真解無共同之處，這樣的差分格式稱為不穩(wěn)定差分格式；反之，如果誤差不在增大或逐步減小，這樣的差分格式稱為穩(wěn)定的差分格式。常用馮.牛曼（Von.Neumann）方法分析差分格式的穩(wěn)定性。*56第五節(jié)差分方程的穩(wěn)定性分析用有限差分方法求解微分方程時的實際誤差為=截斷誤差十舍入誤差

設(shè)時刻的誤差為，時刻的誤差為

差分格式是穩(wěn)定的;

差分格式是不穩(wěn)定的。若*57假設(shè)

時段的誤差為

，差分方程計算解為，則

時刻的真解為

誤差傳播方程第五節(jié)差分方程的穩(wěn)定性分析代真解入上式得:穩(wěn)定性分析:將誤差展成fourier級數(shù)，取其中一項為，并將其寫成：式中：α—某一常數(shù)，可以是實數(shù)，也可以是復(fù)數(shù)

和分別為n+1和n時步的誤差分量，則稱令誤差放大因子誤差傳播方程*60若｜ζ｜≤1，則誤差只能隨時間步的增加而減少，從而保證了差分方程的穩(wěn)定性。稱｜ζ｜≤1為差分方程的穩(wěn)定條件。VonNewmann判別準(zhǔn)則第五節(jié)差分方程的穩(wěn)定性分析*61（1）隱式差分格式的穩(wěn)定性分析誤差傳播方程第五節(jié)差分方程的穩(wěn)定性分析*62由于1所以：｜ζ｜≤1恒成立，由VonNeumann判別準(zhǔn)則可知，隱式差分格式是無條件穩(wěn)定的。第五節(jié)差分方程的穩(wěn)定性分析*63（2）顯示格式的穩(wěn)定條件條件穩(wěn)定（P78）(3)CN格式的穩(wěn)定條件無條件穩(wěn)定第五節(jié)差分方程的穩(wěn)定性分析*64課后編程序練習(xí)一維不穩(wěn)定滲流的顯式求解要求：40個網(wǎng)格，等距空間步長0.2m，不等距時間步長。編制不同初始和邊界條件下壓力分布計算程序。擴展練習(xí)*65第六節(jié)線性代數(shù)方程組的求解方法直接解法迭代解法高斯消元法Lu分解法三對角追趕法D4高斯消元法Jacobi迭代法Gauss-siedel迭代法松弛法(PSOR,LSOR)*66（1）高斯消元法一、直接解法第一步:將方程(1)乘上-1加到方程(3)上,消去(3)中的未知數(shù)x1，可得到*67（1）高斯消元法一、直接解法第二步:將方程(2)加到方程(4)上,消去(4)中的未知數(shù)x2，可得到與原方程組等價的三角形方程組*68（1）高斯消元法一、直接解法a1,1x1+a1,2x2+……+a1,nxn=b1a2,1x1+a2,2x2+……+a2,nxn=b2 ……an,1x1+an,2x2+……+an,nxn=bn下面是n元線性方程組的一般形式：我們可以把它表示為增廣矩陣的形式：a1,1 a1,2 ……

a1,n b1a2,1 a2,2 ……

a2,n b2 ……an,1 an,2 …… an,n bn*69消元過程a1,1(1)

a1,2(1)

……

a1,n(1)

b1(1)a2,1(1)

a2,2(1)

……

a2,n(1)

b2(1) ……an,1(1)

an,2(1)

…… an,n(1)

bn(1)注：用上標(biāo)(k)表示第k次消元前的狀態(tài)第1次消元，第i行的乘數(shù)：（i=2,3,…,n）a1,1(1)

a1,2(1)

……

a1,n(1)

b1(1)

a2,2(2)

……

a2,n(2)

b2(2) ……

an,2(2)

…… an,n(2)

bn(2)得到新的增廣矩陣：ai,j(2)=ai,j(1)-mi,1a1,j(1)bi(2)=bi(1)-mi,1b1(1)（i,j=2,3,…,n）*70第k次消元，第i行的乘數(shù)：（i=k+1,k+2,…,n）消元過程a1,1(1)

a1,2(1)

…………

a1,n(1)

b1(1)

a2,2(2)

…………

a2,n(2)

b2(2) ………… ak,k(k)

…… ak,n(k)

bk(k)

…… an,k(k)

…… an,n(k)

bn(k)第k次消元前的增廣矩陣：第k步消元的主行第k步消元的主元素*71回代過程a1,1(1)

a1,2(1)

……

a1,n(1)

b1(1)

a2,2(2)

……

a2,n(2)

b2(2) …… an,n(n)

bn(n)最后得到的增廣矩陣：最終結(jié)果的計算：高斯消元法是一種比較簡單、適用范圍較廣的有效算法，但在實際應(yīng)用中，我們往往需要具體問題具體分析，對這樣的標(biāo)準(zhǔn)算法進行改進，才能滿足我們的需要。*72由于不能夠保證ai,k(k)是ak,k(k)的倍數(shù)，要想消元，必須使兩行分別乘以一個乘數(shù)。

方程較多時，系數(shù)有可能越來越大，到一定程度有可能導(dǎo)致系數(shù)越界，因此要隨時對各行化簡，即把這一行中所有元素除以這些元素的最大公約數(shù)。

但是，無論如何，這也保證不了不會發(fā)生越界，因此這種算法一般適用于系數(shù)、未知數(shù)范圍較小，未知數(shù)個數(shù)較少的方程。方法分析*73一、直接解法（2）LU（Crout）分解法根本上與高斯消元是一樣的，但是它在分解過程中不會修改到方程右端的向量。且L和U可以保存起來，隨時調(diào)用。可以證明，當(dāng)A的各階順序主子式均不為零時，Doolittle分解可以實現(xiàn)并且唯一。即矩陣A有分解：A=LU，其中L為下三角陣，U為上三角陣。解線性方程組Ax=b可以等價于解兩個三角方程。*74一、直接解法（2）LU分解法*75一、直接解法（2）LU分解法*76一、直接解法（2）LU分解法*77一、直接解法（2）LU分解法*78一、直接解法（2）LU分解法*79（3）三對角追趕法當(dāng)i=1時,bi=0;i=N時,dN=0.一、直接解法一維問題的追趕法求解*80三對角追趕法對應(yīng)矩陣c1d1b2c2d2b3c3d3

bN-1cN-1dN-1

bNcNx1x2x3

xN-1xNg1g2g3

gN-1gN=一、直接解法*81設(shè)代入第二個方程組設(shè)*82設(shè)則可得:按此方法推導(dǎo)消去,可得這幾個系數(shù)的通式:*83方程組的第N-1個方程消元并整理得:設(shè)代上式入第N個方程:*84同樣設(shè):即可得:上述消元過程稱為追趕法的追過程回代xN逐此回代,可得xi的計算公式為:按此公式可得到全部未知量的值.上述消元過程稱為追趕法的趕過程*85一維問題的LU分解法求解*86一維問題的LU分解法求解*87一維問題的LU分解法求解*88其計算工作量為5n-4次乘除法。一維問題的LU分解法求解*89152131516314196174920718231210242111822（4）D4解法一、直接解法152131516314196174920718231210242111822不允許相鄰的網(wǎng)格塊被連續(xù)的置入系數(shù)矩陣中（P66）。*90D4方法的矩陣結(jié)構(gòu)A1A2A3A4一、直接解法*91D4線性代數(shù)方程組的求解AX=BA1A2A3A4X1X2=B1B2一、直接解法*92（1）簡單迭代法N階線性代數(shù)方程組a11x1+a12x2+a13x3++a1NxN=b1a21x1+a22x2+a23x3++a2NxN=b2

aN1x1+aN2x2+aN3x3++aNNxN=bN設(shè)aii不為零,則可將上式變?yōu)榈男问?二、迭代法*93簡單迭代法第六節(jié)線性代數(shù)方程組的求解方法*94給定初始值:代入上式右端可得第一次迭代值:

二、迭代法*95（1）簡單迭代法*96將第一次的值代入迭代式可得:依次下去,可得第K次迭代的值:（1）簡單迭代法*97（1）簡單迭代法-迭代格式*98（1）簡單迭代法-迭代格式*99（2）Gauss-Deidel迭代法*100設(shè)初始值代入上面第一個方程:求出x1(1),用第二個方程求x2(1)時,不再象簡單迭代法那樣單純的將代入,而是用新算出的x1(1)代替x1(0),即（2）Gauss-Deidel迭代法*101（2）Gauss-Deidel迭代法-迭代格式*102（2）Gauss-Deidel迭代法-迭代格式*103（3）超松弛迭代法(SOR)