地下水?dāng)?shù)值模擬03

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

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

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

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

5、()(dxfdxdxfdxdxdfxxfxxf差分公式對(duì)比名稱名稱公式公式截?cái)嗾`差截?cái)嗾`差一一階階導(dǎo)導(dǎo)數(shù)數(shù)前差前差后差后差中心差中心差二階導(dǎo)數(shù)二階導(dǎo)數(shù)xxfxxfdxdf)()(xxxfxfdxdf)()(xxxfxxfdxdf2)()()(2xO )( xO )( xO 222)()()(2)(xxxfxfxxfdxfd)(2xO 思考與練習(xí) 高階導(dǎo)數(shù)的有限差分形式高階導(dǎo)數(shù)的有限差分形式33dxfd?yxf2?44xf? dxxdfdxxxdfdxxxdfxxxfxxfxxfdxddxfddxddxfd21 2222233 332222222222122221xxxfxxfxxfxxfxxf

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

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

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

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

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

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

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

13、tT定義定義kikikikihhhh111)21 (截?cái)嗾`差為：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對(duì)i=0,1,n循環(huán)執(zhí)行xiHhi00t=t+t對(duì)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)(整理得：整理得：截?cái)嗾`差為：O(x2)+O(t) 2*)( xtT定義定義kikikikihhhh11111)21 (hik+1，。hik+1（i=1,i=2,i=l-1)差分方程代入形成由（l-1)個(gè)方程組成的線性方程組，聯(lián)立求解。 1122101112121112

15、121212knknknkkknknkkhhhhhhhhhh開始輸入初始參數(shù)t=0對(duì)i=0,1,n循環(huán)執(zhí)行xiHhi00t=t+t對(duì)i=1,n-1循環(huán)線性方程組調(diào)用追趕法程序求解三對(duì)角線方程組輸出結(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時(shí)刻k+1時(shí)刻2t(i,k)k+1/2時(shí)刻虛擬時(shí)刻變量虛擬時(shí)刻變量怎么?。吭趺慈?？tHxHT*22截?cái)嗾`差為：O(x2)+O (t2 ) 21*kitHkikikixHxHTxHT2212221222121中心差分中心差分thhxhhhxhhhTkikikikikikikiki1*211211111)(2)(22整理得：整理得：kikikikikikihhhhhh1111111)22()22(2*)( xtT定義定義hik+1，thhxhhhxhhhkikikikikikikiki1211211111T)(21)(2）（ 小結(jié)顯式差分隱式差分中心式差分差分方程公式向前差分向前差分向后差分向后差分中心