地下水數(shù)值模擬03

1、第三章第三章 有限差分法有限差分法主要內(nèi)容主要內(nèi)容 有限差分法的基本原理有限差分法的基本原理 幾種主要的差分格式幾種主要的差分格式 二維滲流問題的差分方程二維滲流問題的差分方程 一般差分方程組的解法一般差分方程組的解法第一節(jié)第一節(jié)有限差分法的有限差分法的基本原理基本原理一、有限差分法的基本思想一、有限差分法的基本思想 1、基本原理、基本原理 從物理現(xiàn)象引出相應從物理現(xiàn)象引出相應微分方程微分方程（方程（方程+邊界條件）；邊界條件）； 用差分網(wǎng)格用差分網(wǎng)格離散離散求解域；求解域； 用差分公式將基本方程轉(zhuǎn)化為用差分公式將基本方程轉(zhuǎn)化為差分方程差分方程（代數(shù)方程）；（代數(shù)方程）； 用差分方程的解作為微

2、分方程的近似解。用差分方程的解作為微分方程的近似解。一、有限差分法的基本思想一、有限差分法的基本思想 2、求解步驟、求解步驟選取網(wǎng)格；選取網(wǎng)格；對微分方程及定解條件選擇差分近似，列出對微分方程及定解條件選擇差分近似，列出差分格式；差分格式；求解差分格式；求解差分格式；討論差分格式解對于微分方程解的收斂性及討論差分格式解對于微分方程解的收斂性及誤差估計。誤差估計。差分網(wǎng)格差分網(wǎng)格xyx，y空間步長空間步長t時間步長時間步長結(jié)點格點一、有限差分法的基本思想一、有限差分法的基本思想 3、優(yōu)缺點、優(yōu)缺點優(yōu)優(yōu) 點點缺缺 點點 1. 1.簡單問題的數(shù)學表達簡單問題的數(shù)學表達式和計算的執(zhí)行過程式和計算的執(zhí)行

3、過程比較直觀、易懂；比較直觀、易懂； 2.2.算法效率比較高，易算法效率比較高，易編程；編程；1 1、對自然邊界處理的靈、對自然邊界處理的靈活性較差?；钚暂^差。2 2、對溶質(zhì)運移等問題，、對溶質(zhì)運移等問題，精度受限精度受限二、導數(shù)的有限差分近似表示二、導數(shù)的有限差分近似表示 1、差分的概念、差分的概念xxfxxfxydxdyxx)()(limlim00設有x的解析函數(shù)y(x)，函數(shù)y對x的導數(shù)為： 是函數(shù)對自變量的導數(shù)，又稱微商。是函數(shù)對自變量的導數(shù)，又稱微商。dxdy 、 分別稱為函數(shù)及其自變量的差分分別稱為函數(shù)及其自變量的差分yxdy、dx分別是函數(shù)及自變量的微分分別是函數(shù)及自變量的微分,

4、xy 為函數(shù)對自變量的差商。為函數(shù)對自變量的差商。由導數(shù)（微商）和差商的定義可知，由導數(shù)（微商）和差商的定義可知，當自變量的差分（增量）趨近于零時，當自變量的差分（增量）趨近于零時，就可以由差商得到導數(shù)。因此在數(shù)值就可以由差商得到導數(shù)。因此在數(shù)值計算中常用差商近似代替導數(shù)。計算中常用差商近似代替導數(shù)。二、導數(shù)的有限差分近似表示二、導數(shù)的有限差分近似表示 2、導數(shù)的有限差分形式、導數(shù)的有限差分形式用泰勒級數(shù)展開可以推導出導數(shù)的有限差分形式。用泰勒級數(shù)展開可以推導出導數(shù)的有限差分形式。333222! 3)(! 2)()()(dxfdxdxfdxdxdfxxfxxf333222! 3)(! 2)()

5、()(dxfdxdxfdxdxdfxxfxxf差分公式對比名稱名稱公式公式截斷誤差截斷誤差一一階階導導數(shù)數(shù)前差前差后差后差中心差中心差二階導數(shù)二階導數(shù)xxfxxfdxdf)()(xxxfxfdxdf)()(xxxfxxfdxdf2)()()(2xO )( xO )( xO 222)()()(2)(xxxfxfxxfdxfd)(2xO 思考與練習 高階導數(shù)的有限差分形式高階導數(shù)的有限差分形式33dxfd?yxf2?44xf? dxxdfdxxxdfdxxxdfxxxfxxfxxfdxddxfddxddxfd21 2222233 332222222222122221xxxfxxfxxfxxfxxf

6、xxfxxfxxfxxxxfxxfxxxfxfxxfxxfx思考與練習 高階導數(shù)的有限差分形式高階導數(shù)的有限差分形式y(tǒng)xf2?867527865310241 2221 2ffffhhffhffhhyfyfyxf015234678910111211931042031231102109220223221220444461 22221 2fffffhhfffhfffhfffhhxfxfxfxf三、簡單水文地質(zhì)模型的有限差分方程組三、簡單水文地質(zhì)模型的有限差分方程組 1、水文地質(zhì)模型、水文地質(zhì)模型u 以一維河間地塊承壓含水層中的水流問題為例。以一維河間地塊承壓含水層中的水流問題為例。u 考慮一個以通過

7、考慮一個以通過x0和和xL處的長且直的河流為界的承壓含水層，該含水處的長且直的河流為界的承壓含水層，該含水層均質(zhì)各向同性，頂?shù)装逅?，上覆弱透水層，垂向補給強度為層均質(zhì)各向同性，頂?shù)装逅?，上覆弱透水層，垂向補給強度為W，兩河流，兩河流邊界的水位分別邊界的水位分別0、 L為且不隨時間變化。為且不隨時間變化。u 試研究含水層的水頭分布。試研究含水層的水頭分布。 三、簡單水文地質(zhì)模型的有限差分方程組三、簡單水文地質(zhì)模型的有限差分方程組 2、數(shù)學模型、數(shù)學模型LLxxxHxHLxWxHT)()()0( 00022212222CxCxTWHCxTWxHTWxHLLCLTWLx0122時，當020Cx時

8、，當02012CLTWLCL 00222xTWLLxTWxHL三、簡單水文地質(zhì)模型的有限差分方程組三、簡單水文地質(zhì)模型的有限差分方程組 3、有限差分方程、有限差分方程LLxxxHxHLxWxHT)()()0( 00022（1）網(wǎng)格剖分）網(wǎng)格剖分 沿河流的方向取單寬0,L作為研究區(qū)域，將L等分為n份，空間步長 對每個結(jié)點進行編號，結(jié)點編號由左向右依次為0，1，2，i，n。共有n+1個結(jié)點，其中2個邊界結(jié)點，n-1個內(nèi)結(jié)點。nLx 三、簡單水文地質(zhì)模型的有限差分方程組三、簡單水文地質(zhì)模型的有限差分方程組 3、有限差分方程、有限差分方程（2）建立有限差分方程）建立有限差分方程先任取一結(jié)點i進行分析。

9、 022WxHT0)(2211WxhhhTiii2112xTWhhhiii移項整理，得：移項整理，得： 對于所有內(nèi)結(jié)點內(nèi)結(jié)點1、2、n-1，建立結(jié)點相應的差分方程組 LnnnnnTxWhhTxWhhhTxWhhhTxWhh2122123232102212 2 2 2n-1個線性代數(shù)方程，未知量共n-1個故方程可解。 三、簡單水文地質(zhì)模型的有限差分方程組三、簡單水文地質(zhì)模型的有限差分方程組 3、有限差分方程、有限差分方程（3）求解）求解LnnnnnTxWhhTxWhhhTxWhhhTxWhh2122123232102212 2 2 2LnnTxWTxWTxWTxWhhhh/21121121122

10、22021221將方程表示成矩陣形式，則有： 線性代數(shù)方程組！線性代數(shù)方程組！ 系數(shù)矩陣中的元素僅位于三條對角線上，系數(shù)矩陣對稱且系數(shù)矩陣中的元素僅位于三條對角線上，系數(shù)矩陣對稱且正定，故稱為正定，故稱為三對角線性代數(shù)方程組。三對角線性代數(shù)方程組。 最有效的求最有效的求解方法解方法追趕法追趕法第二節(jié)第二節(jié)幾種主要的差分格式幾種主要的差分格式水文地質(zhì)模型描述水文地質(zhì)模型描述u 以一維河間地塊承壓含水層中的水流問題為例。u 含水層均質(zhì)各向同性，不考慮垂向補給，兩河流邊界的水位隨時間變化，分別0(t)、 L(t) 。u 試研究含水層的水頭分布。 )0( )(),()0( )(),()0( )(),(

11、)0 ,0( 0000*22TtttxHTtttxHLxxHtxHTtLxtHxHTLLxxt該問題屬于一維承壓非穩(wěn)定流的定解問題。求解該問題，需要對空間該問題屬于一維承壓非穩(wěn)定流的定解問題。求解該問題，需要對空間和時間進行離散，形成的差分網(wǎng)格稱為時空網(wǎng)格。和時間進行離散，形成的差分網(wǎng)格稱為時空網(wǎng)格。 網(wǎng)格剖分網(wǎng)格剖分以等距剖分為例以等距剖分為例將研究區(qū)域?qū)⒀芯繀^(qū)域0，L用直線等分為用直線等分為n份，把時間段份，把時間段0,T用直線等分成用直線等分成m份份以以 表示結(jié)點（表示結(jié)點（i,k）處的水頭）處的水頭 kiH)0( )(),()0( )(),()0( )(),()0 ,0( 0000*2

12、2TtttxHTtttxHLxxHtxHTtLxtHxHTLLxxt 導數(shù)可以利用一階、二階導數(shù)的差商代替，由于一階導數(shù)可以有三種差商表示，因此分別對水頭關于時間的導數(shù)項分別運用前差、后差、中心差將得到三種差分格式。 顯式有限差分顯式有限差分 前差前差 隱式有限差分隱式有限差分 后差后差 中心式有限差分中心式有限差分 中心差中心差 一、一維顯式有限差分格式一、一維顯式有限差分格式tHxHT*22(i,k)(i-1,k)(i+1,k)(i,k+1)211)(2xhhhTkikikithhkiki1*向前差分向前差分=kikikikikihhhhhxtT1112*2)(整理得：整理得：2*)( x

13、tT定義定義kikikikihhhh111)21 (截斷誤差為：O(x2)+O(t) ktkhikktk+1hik+1，。 h11， h21 hn-11h01hl1,差分公（取k=1， i=1，i=2,i=n-1，t2t3 ，t4 ，t0h00， h10，hl0；差分方程，k=0，i=1，i=2,i=n-1,t1h11，hn-10。開始輸入初始參數(shù)t=0對i=0,1,n循環(huán)執(zhí)行xiHhi00t=t+t對i=1,n-1循環(huán)執(zhí)行kikikikihhhh111)21 ( th010 thLl1輸出結(jié)果：t，1ih10iihh tTsum?結(jié)束NoYes hik+1Hik+1?顯式差分格式收顯式差分格

14、式收斂和穩(wěn)定的條件斂和穩(wěn)定的條件2102)(Txt ，x, t：x, t：二、一維隱式有限差分格式二、一維隱式有限差分格式tHxHT*22(i,k+1)(i-1,k+1)(i+1,k+1)(i,k)211111)(2xhhhTkikikithhkiki1*向后差分向后差分=kikikikikihhhhhxtT1111112*2)(整理得：整理得：截斷誤差為：O(x2)+O(t) 2*)( xtT定義定義kikikikihhhh11111)21 (hik+1，。hik+1（i=1,i=2,i=l-1)差分方程代入形成由（l-1)個方程組成的線性方程組，聯(lián)立求解。 1122101112121112

15、121212knknknkkknknkkhhhhhhhhhh開始輸入初始參數(shù)t=0對i=0,1,n循環(huán)執(zhí)行xiHhi00t=t+t對i=1,n-1循環(huán)線性方程組調(diào)用追趕法程序求解三對角線方程組輸出結(jié)果：t，1ih10iihh tTsum?結(jié)束NoYeskikikikihhhh11111)21 ( 三、一維中心式有限差分格式三、一維中心式有限差分格式tk與與tk+1tk+t/2（i,k+1/2)(i,k+1/2)(i-1,k+1/2)(i+1,k+1/2)(i,k+1)k時刻k+1時刻2t(i,k)k+1/2時刻虛擬時刻變量虛擬時刻變量怎么??？怎么取？tHxHT*22截斷誤差為：O(x2)+O (t2 ) 21*kitHkikikixHxHTxHT2212221222121中心差分中心差分thhxhhhxhhhTkikikikikikikiki1*211211111)(2)(22整理得：整理得：kikikikikikihhhhhh1111111)22()22(2*)( xtT定義定義hik+1，thhxhhhxhhhkikikikikikikiki1211211111T)(21)(2）（ 小結(jié)顯式差分隱式差分中心式差分差分方程公式向前差分向前差分向后差分向后差分中心