地下水流數值模擬方法第六章

水動力彌散方程的數值解法中國地質大學環境學院2019春地下水流數值模擬方法第六章水動力彌散方程的數值解法中國地質12一、有限差分法基本思想：

將空間劃分成若干網格，把時間分成若干小時段，每一個網格中心點（格點）處的未知變量（H或C）視為該網格平均值。利用差商代替微商。基本步驟：

（1）剖分（2）建立差分方程組（3）求解2一、有限差分法基本思想：23一、有限差分法-導數的有限差分近似圖中，去x軸上任意一點i，其坐標為

在改點左右相聚為處分別取（i-1）和（i+1），其坐標分別為和以i為中心，泰勒展開C(x)3一、有限差分法-導數的有限差分近似圖中，去x軸上任意一點i34一、有限差分法-導數的有限差分近似整理并略去余項（6-1）-（6-2），再除以略去余項（6-1）-（6-2），再除以略去余項：分別一階導數的向前差分、向后差分及中心差分二階中心差分4一、有限差分法-導數的有限差分近似整理并略去余項分別一階導45一、有限差分法-一維水動力彌散的差分解法一維水動力彌散方程

（6-7）（1）顯格式式中僅有一個未知數，解得

式（6-7）中的對流項取中心差分5一、有限差分法-一維水動力彌散的差分解法一維水動力彌散方程56

可以證明，穩定性準則要求

即

即。（1）（2）。由條件（2），格距要求很小；由（2）可知，鑒于較小，導致不能太大，將（2）代入（1）式中，得到

6

可以證明，穩定性準則要求

即67若對流項改為向后差分

解得

穩定性要求不難看出，穩定性限制比對流項取中心差分有所改善。7若對流項改為向后差分解得78（2）隱式格式

整理后得到

隱格式是無條件穩定的。8（2）隱式格式89（2）Grank-Nicolson格式

整理后得到

取隱式和顯示的平均，即Grank-Nicolson格式9（2）Grank-Nicolson格式910一、有限差分法-二維水動力彌散的差分解法二維水動力彌散方程

（4-56）（1）顯格式式中僅有一個未知數式（4-56）中的對流項取中心差分10一、有限差分法-二維水動力彌散的差分解法二維水動力彌散方1011一、有限差分法-二維水動力彌散的差分解法化簡后，有涉及以（i，j）為中心的5個網格點在tn時刻的已知濃度11一、有限差分法-二維水動力彌散的差分解法化簡后，有涉及以1112一、有限差分法-二維水動力彌散的差分解法（2）隱格式式（4-56）中右端的對流項取中心差分，右端個C的時階均取n+1水平12一、有限差分法-二維水動力彌散的差分解法（2）隱格式式（1213一、有限差分法-二維水動力彌散的差分解法（2）隱格式整理后收斂且無條件穩定涉及以（i，j）為中心的5個網格點在tn+1時刻的未知濃度13一、有限差分法-二維水動力彌散的差分解法（2）隱格式整理1314一、有限差分法-二維水動力彌散的差分解法（3）Grank-Nicolson格式將隱式格式的兩式相加除以214一、有限差分法-二維水動力彌散的差分解法（3）Grank1415一、有限差分法-二維水動力彌散的差分解法（3）Grank-Nicolson格式整理得15一、有限差分法-二維水動力彌散的差分解法（3）Grank1516一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱式法（ADI法）分兩次對三對角矩陣求逆，將一個Δt分成兩個Δt/2計算第一個半時間步，對x方向的偏導數采用隱式差分，對y方向的偏導數采用顯示差分。16一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱1617一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱式法（ADI法）整理得第二個半時間步，對y方向的偏導數采用隱式差分，對x方向的偏導數采用顯示差分。17一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱1718一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱式法（ADI法）整理得收斂且無條件穩定18一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱1819一、有限差分法-求解差分方程的計算機程序舉例算例

見教材P81-8519一、有限差分法-求解差分方程的計算機程序舉例算例1920二、有限單元法-伽遼金有限單元法伽遼金法

對均勻多孔介質，一維穩定流場中二維水動力彌散方程

取x軸方向與地下水流向相同，記伽遼金法是尋找一個級數形式的試探函數作微分方程的近似解，并使其滿足邊界條件

（6-26）20二、有限單元法-伽遼金有限單元法伽遼金法（6-26）2021二、有限單元法-伽遼金有限單元法上式一般不會滿足方程，因為僅是近似解，得到一個剩余誤差函數R(x，y)，在平均意義下迫使誤差為0

我們加以權，令剩余的加權積分為0，W是一組權函數。加權剩余法（6-29）21二、有限單元法-伽遼金有限單元法上式一般不會滿足方程，因2122二、有限單元法-伽遼金有限單元法在整個研究區D上，基函數NL(x,y)分段定義(1)函數區域連續（2）一階導數不連續（3）二階導不容易確定故采用分步積分的方法使二階導數降階22二、有限單元法-伽遼金有限單元法在整個研究區D上，基函數2223二、有限單元法-伽遼金有限單元法對一維積分整理后分步積分推廣到二維的情況（6-32）代入（6-29）23二、有限單元法-伽遼金有限單元法對一維積分（6-32）代2324二、有限單元法-伽遼金有限單元法有限單元剖分與基函數（1）三角形單元見《地下水流動問題數值模擬》（2）矩形單元將區域記為D，邊界記為B。要求各單元均質、等厚，即T、μ為常數。結點（內結點、邊界結點）

24二、有限單元法-伽遼金有限單元法有限單元剖分與基函數2425二、有限單元法-伽遼金有限單元法構造單元函數NL基函數取為“雙線性插值”

基函數NL（x，y）形狀如同一頂高等于1、有4條直線斜邊和4條下凹型曲邊的尖頂斗笠25二、有限單元法-伽遼金有限單元法構造單元函數NL基函數N2526二、有限單元法-伽遼金有限單元法子區DL以結點L為其共同結點的所有矩形單元（<=4）基函數NL僅在子區上不為0，在非DL上均為0,屬于子區DL單元e的單元基函數NL，用NeL表示（>0）26二、有限單元法-伽遼金有限單元法子區DL2627二、有限單元法-伽遼金有限單元法典型矩形單元該單元中心位于坐標原點處，且4條邊分別平行x，y軸,結點從左下角開始按逆時針方向編號i,j,m,n。沿x方向長2a，y方向寬2b。27二、有限單元法-伽遼金有限單元法典型矩形單元2728二、有限單元法-伽遼金有限單元法按上述要求所構造的NeL形式為對典型單元，令

則28二、有限單元法-伽遼金有限單元法按上述要求所構造的NeL2829二、有限單元法-伽遼金有限單元法其中即（6-34）29二、有限單元法-伽遼金有限單元法其中即（6-34）2930二、有限單元法-伽遼金有限單元法Nei(x,y)在結點i處為1，其它3處為0根據上式，有平行x或y方向上，Nei(x,y)線性變化故單元基函數為雙線性插值函數若結點坐標用表示結點濃度用表示30二、有限單元法-伽遼金有限單元法Nei(x,y)在結點i3031二、有限單元法-伽遼金有限單元法結點濃度CL(t)和單元基函數NeL(x,y)來確定單元內的近似解，即對于區域D，結點L的基函數在子區DL內各單元分塊確定。稱為矩陣單元e上的基函數（6-36）31二、有限單元法-伽遼金有限單元法結點濃度CL(t)和單元3132二、有限單元法-伽遼金有限單元法伽遼金有限單元法

左端積分范圍為D，由于在上NL=0,改寫在子區DL的積分，對矩形單元逐個積分、求和。設子區DL上有neL個單元，有（6-39）32二、有限單元法-伽遼金有限單元法伽遼金有限單元法（6-33233二、有限單元法-伽遼金有限單元法由（6-36）得由（6-34）得再由（6-36）得（6-40）（6-42）33二、有限單元法-伽遼金有限單元法由（6-36）得由（6-3334二、有限單元法-伽遼金有限單元法將（6-40）（6-42）代入（6-39）得34二、有限單元法-伽遼金有限單元法將（6-40）（6-43435二、有限單元法-伽遼金有限單元法（6-44）35二、有限單元法-伽遼金有限單元法（6-44）3536二、有限單元法-伽遼金有限單元法用矩陣形式表示一階微分方程組將一階導數用差分形式表示按矩陣符號有若﹛C﹜在新時間水平t+Δt上取值，則36二、有限單元法-伽遼金有限單元法用矩陣形式表示一階微分方3637二、有限單元法-伽遼金有限單元法（1）給定初始條件和邊界條件，可求（6-46）各個系數矩陣，于是可解得第一個末時刻各結點的濃度；（2）以此濃度為已知濃度，計算新的系數矩陣；（3）再次求解，得到第二個時間步長末時刻結點濃度（6-46）濃度可在之間任意位置近似表示37二、有限單元法-伽遼金有限單元法（1）給定初始條件和邊界3738二、有限單元法-伽遼金有限單元法按單元順序計算系數矩陣每個矩陣單元涉及4個結點i，j，m，n，涉及4個方程單元e的作用，可用4行與4列的矩陣來表示，稱為單元矩陣38二、有限單元法-伽遼金有限單元法按單元順序計算系數矩陣3839二、有限單元法-伽遼金有限單元法算例：單元彌散矩陣與整體彌散矩陣的關系設研究區劃分為8個單元，16個結點，四周為隔水邊界編號如下：39二、有限單元法-伽遼金有限單元法算例：單元彌散矩陣與整體3940二、有限單元法-伽遼金有限單元法算例：單元彌散矩陣與整體彌散矩陣的關系計算每個單元彌散矩陣，按雙下標已知，放置單元彌散矩陣并疊加，得到總體彌散矩陣。40二、有限單元法-伽遼金有限單元法算例：單元彌散矩陣與整體4041二、有限單元法-伽遼金有限單元法從（6-44）知將坐標原點移到矩形單元中心有41二、有限單元法-伽遼金有限單元法從（6-44）知將坐標原4142二、有限單元法-伽遼金有限單元法通過雙下標一致的元素求和，得到總體彌散矩陣諸元素42二、有限單元法-伽遼金有限單元法通過雙下標一致的元素求和4243二、有限單元法-伽遼金有限單元法積分求解可用高斯求積方法，詳見P97-10543二、有限單元法-伽遼金有限單元法積分求解可用高斯求積方法4344二、有限單元法-伽遼金有限單元法44二、有限單元法-伽遼金有限單元法4445二、有限單元法-伽遼金有限單元法算例：二維水動力彌散的伽遼金法計算機程序

數學模型寫成剖分成9個單元共20個結點，空間步長取Δx=Δy=10m，時間步長Δt=5d，滲透流速u=0.01m/d具體代碼見P106-11145二、有限單元法-伽遼金有限單元法算例：二維水動力彌散的伽4546三、過程與數值彌散過量現象

如一維流動水動力彌散中單位濃度梯度函數注入的情況。

對比解析解與數值解C-x曲線，在濃度封面附近數值計算的濃度超過最大濃度值1和小于最小濃度0，失去物理意義。數值彌散

若純對流方程采用有限差分逼近，其結果呈現似有彌散的存在。46三、過程與數值彌散過量現象數值彌散4647三、過程與數值彌散47三、過程與數值彌散4748三、過程與數值彌散過量現象-解析

由于時間步長與空間尺度匹配不佳，因而含水層在數值上不能“吸收”住污染物所引起。解決辦法：

將時間增量選擇為幾何級數的子項；

做試驗校正時間步長和差分網格；48三、過程與數值彌散過量現象-解析4849三、過程與數值彌散數值彌散-解析

由截斷誤差引起。

對流項一階空間導數采用向后差分逼近，濃度Ci-1泰勒展開變形得

向后差分的截斷誤差主部分為相當于彌散系數的彌散項49三、過程與數值彌散數值彌散-解析相當于彌散系數的彌散項4950三、過程與數值彌散數值彌散-解析

由故純對流方程有限差分近似表示后，結果如同對流-彌散方程的結果，記用數值彌散系數與物理彌散系數的比值來評價數值彌散對物理彌散的干涉程度此條件下的數值彌散系數50三、過程與數值彌散數值彌散-解析此條件下的數值彌散系數5051三、過程與數值彌散數值彌散-解析

定義為Peclet數，時，彌散占優；

時，對流占優。速度u與物理彌散系數DL比越大，越大。對于小數，數值解與精確解非常接近，對于大數，數值解過渡帶變寬，濃度鋒面變平緩，在鋒面的前后還出現解的振動—過量現象。控制Δx的大小從而減少數值彌散51三、過程與數值彌散數值彌散-解析控制Δx的大小從而減少數5152三、過程與數值彌散對純對流方程

52三、過程與數值彌散對純對流方程5253三、過程與數值彌散對純對流方程

結果也造成數值彌散此時，總的數值彌散是兩者之和53三、過程與數值彌散對純對流方程結果也造成數值彌散此時，總5354四、對流占主導地位時的數值方法上游加權法

對流項的有限差分化過程中乘上一個適當的權因子，則波動和過量得到改善或消除一維水動力彌散中對流項取中心隱式差分若相鄰格點i-1或（和）i+1的濃度越大，則格點i的濃度也越大，反映到差分方程上，系數Gi-1和Gi+1不得與Gi同號，上式要求

不能滿足時，對流占優勢54四、對流占主導地位時的數值方法上游加權法對流項的有限差分5455四、對流占主導地位時的數值方法上游加權法

假定DL和u為常數，取等格距差分網格。一維對流-彌散方程在積分，有對流-彌散通量用有限差分近似表示

55四、對流占主導地位時的數值方法上游加權法假定DL和u為常5556四、對流占主導地位時的數值方法上游加權法

ω為權因子；當ω=1/2時，對流項相當于中心差分；當ω=1時，對流項相當于取向后差分，

當ω=0時，對流項相當于向前差分。速度不同，

ω取值不同56四、對流占主導地位時的數值方法上游加權法5657四、對流占主導地位時的數值方法上游加權法

上游濃度權因子大于1/2。物理意義：加強上游格點濃度值在對流項中的作用可通過濃度的線性（或二次）插值具體確定在Δt期間C在格邊上隨時間的變化，從而確定ω值。57四、對流占主導地位時的數值方法上游加權法上游濃度權因子大5758四、對流占主導地位時的數值方法

仍用Gi-1、Gi和Gi+1分別表示上式左端第一、二、三項括號內系數。則采用上游加權法有限差分格式，可消除過量和波動，但截斷誤差增大，數值彌散反而增大。中心差分格式，截斷誤差為二階上游加權有限差分截斷誤差為一階58四、對流占主導地位時的數值方法仍用Gi-1、Gi5859四、對流占主導地位時的數值方法

上游權因子的物理實質：對流項濃度C在某斷面，例如斷面的中心差分意味著對流項的濃度在Δt時段內均以i與i-1兩格點的平均濃度，即通過斷面；實際為：Δt時段的初始時刻是該平均濃度值，隨后在對流的作用下，該斷面的濃度顯然會向上游格點的濃度Ci-1(u>0時)的變化。59四、對流占主導地位時的數值方法上游權因子的物理實5960四、對流占主導地位時的數值方法特征值方法

令源匯項為0，展開對流-彌散方程

60四、對流占主導地位時的數值方法特征值方法6061四、對流占主導地位時的數值方法特征值方法

假定流體密度值ρ不變（不可壓縮），含水層不可壓縮，則，（6-75）簡化為基本思路：物質導數表示一個沿軌跡（特征曲線）

運動的研究對象濃度變化率。由彌散作用引起。若存在運會想，還將受到沿軌跡源與匯、化學、生物化學反應、吸附與解析等。對流項用沿特征線的示蹤質點的運動描述

（6-76）61四、對流占主導地位時的數值方法特征值方法（6-76）6162四、對流占主導地位時的數值方法特征值方法

特征點均勻分布研究區，賦一個初始濃度。

含水層中未受污染部分特征點的初始濃度為零，特征點濃度眼軌跡的瞬時濃度變化，借助矩形網格差分格式計算。62四、對流占主導地位時的數值方法特征值方法6263四、對流占主導地位時的數值方法63四、對流占主導地位時的數值方法6364四、對流占主導地位時的數值方法模型精確度取決于網格距大小。可以在同一網格上計算速度場。如果采用較少的網格距，可以改善流場的描述。

坐標為的特征點所在的單元，可用下標i,j確定

精確度受到每一網格內所選擇的最初的特征點數目影響。64四、對流占主導地位時的數值方法模型精確度6465四、對流占主導地位時的數值方法缺點：（1）實際程序冗長，需記錄特征點的蹤跡；

（2）抽水點或出流邊界的特征點需要消除，入滲邊界須有特征點生成，不透水邊界特征點須反射回來；

（3）要通過消去特征點來避免滯留點附近特征點堆積；（4）特征點耗盡的單元，必須產生新的特征點。65四、對流占主導地位時的數值方法缺點：6566四、對流占主導地位時的數值方法動坐標系方法與網格變形方法能消除過量現象，還保持較陡的濃度鋒面。

對式

設坐標變換

得，

取空間步長，則

66四、對流占主導地位時的數值方法動坐標系方法與網格變形方法6667四、對流占主導地位時的數值方法動坐標系方法與網格變形方法

在固定坐標系中取格距為，有

動坐標系中格點運動與固定坐標系格點運動重合。

設固定坐標系中格點xj在tn+1時刻與動坐標系中Xi重合，有

動坐標系下濃度結點的表達式67四、對流占主導地位時的數值方法動坐標系方法與網格變形方法6768四、對流占主導地位時的數值方法動坐標系方法與網格變形方法

式變成

對半無限長砂柱，當DL=0時68四、對流占主導地位時的數值方法動坐標系方法與網格變形方法6869四、對流占主導地位時的數值方法網格變形方法

設tn時刻有限元結點在空間中的分布是對應濃度分布為

69四、對流占主導地位時的數值方法網格變形方法6970四、對流占主導地位時的數值方法網格變形方法

主要步驟：結點位置隨時間變化，基函數及有限元方程的系數都依賴于時間，每推進一個Δt都需要重新計算。二維問題還須進行畸形判斷

70四、對流占主導地位時的數值方法網格變形方法7071四、對流占主導地位時的數值方法結合動坐標的網格變形方法

一維對流-彌散在某個動點在t時刻位置x可寫成

其中，x0是動點的初始位置，有取時間的一階偏導數對流-彌散方程寫成

若控制動點速度dx/dt并讓它接近流速u，則可轉換成彌散為主的的方程式71四、對流占主導地位時的數值方法結合動坐標的網格變形方法若7172四、對流占主導地位時的數值方法若使用伽遼金法，基函數因依賴于結點位置而成為時間的函數。若用表示隨時間變化的基函數，則對任意t，有

故72四、對流占主導地位時的數值方法若使用伽遼金法，基函數因依7273四、對流占主導地位時的數值方法隨機步行法

采用示蹤劑描述污染物運移。只有被污染的特征點在流場中運動。每個特征點都給定一個不變的污染物質量，此質量的和等于排進含水層的總污染物質量。取許多單個特征點軌跡（隨機步行）的平均值，可得到一個有彌散的特征點的分布。要得到一個濃度分布，就要先疊加一個網格并算每個網格單元所含有的特征點數。

73四、對流占主導地位時的數值方法隨機步行法7374四、對流占主導地位時的數值方法隨機步行法-以一維運移為例在x=0處，一個瞬時注入的、質量為Δ

M的理想失蹤及，在t時刻濃度分布對一個固定時間t，看做正態分布74四、對流占主導地位時的數值方法隨機步行法-以一維運移為例7475四、對流占主導地位時的數值方法隨機步行法-以一維運移為例在時間t=0和位置x=0處，放置具有污染物質量為ΔM/N的N個特征點，每一個特征點移動距離x，而在時刻t到達它們各自位置，即Z是一個平均值為0，標準差為1的正態分布隨機變量

所得的頻率分布，通過一個標準化因子，有對有限數量的顆粒，可近似去頂C(x,t)，將空間變量劃分為長度Δx的若干區間，區間中點濃度75四、對流占主導地位時的數值方法隨機步行法-以一維運移為例7576四、對流占主導地位時的數值方法隨機步行法-以一維運移為例若孔隙平均流速是時間和位置的函數，那么時間為0開始的顆粒P在時間間隔Δt內的隨機軌跡xp(t)為

76四、對流占主導地位時的數值方法隨機步行法-以一維運移為例7677四、對流占主導地位時的數值方法隨機步行法-二維流場

首先假定流動方向平行于x軸，得到Z和Z’是正態分布的隨機變量的兩個值。對任意方向的流動，需考慮彌散張量特性。77四、對流占主導地位時的數值方法隨機步行法-二維流場7778四、對流占主導地位時的數值方法隨機步行法-二維流場

首先假定流動方向平行于x軸，得到Z和Z’是正態分布的隨機變量的兩個值。對任意方向的流動，需考慮彌散張量特性。對流運動可確定在水流流動方向及其正交方向的彌散運動78四、對流占主導地位時的數值方法隨機步行法-二維流場對流運7879四、對流占主導地位時的數值方法隨機步行法

79四、對流占主導地位時的數值方法隨機步行法7980四、對流占主導地位時的數值方法隨機步行法

數值結果的質量主要取決于特征點數目N，要考慮網格的選擇。大網格間距會產生極平均的結果，小網格間距會產生十分粗糙的分布。

在穩定流動和一個源或數個源都有固定強度的情況下，可以節省很多特征點。80四、對流占主導地位時的數值方法隨機步行法8081818182四、對流占主導地位時的數值方法隨機步行法

nδ是在時間τ內，有在時間t=0時瞬時注入所產生的分布。所得的特征點分布加權疊加到最終的特征點分布上。用遲滯流速u/R代替平均孔隙流速u，并用遲滯系數去除全部注入質量，可引入線性吸附作用。分配給特征點一個按時間變化的污染物質量82四、對流占主導地位時的數值方法隨機步行法8283四、對流占主導地位時的數值方法隨機步行法-特點

不會出現數值彌散，但在污染帶邊緣的得到滿意的計算結果需要大量的特征點。

完整的，易于變成并能用于每一個水流模型的方法，易推廣到三維模型。

靈敏度較小的參數的變化引起的濃度分布變化會被隨機性濃度的變化覆蓋。

83四、對流占主導地位時的數值方法隨機步行法-特點8384四、對流占主導地位時的數值方法引入人工擴散量方法

對流-彌散方程簡寫成簡寫成腳標注釋的形式引入人工擴散系數并賦權得

（6-107）84四、對流占主導地位時的數值方法引入人工擴散量方法（6-18485四、對流占主導地位時的數值方法引入人工擴散量方法

當θu=0時，導出對稱的正定矩陣，若不校正，格式是不穩定的。在時間步長中點泰勒展開

（6-109）85四、對流占主導地位時的數值方法引入人工擴散量方法（6-18586四、對流占主導地位時的數值方法引入人工擴散量方法

（6-107）對時間求偏導并代入（6-109）得在時，是具有二階截斷誤差的Grank-Nicolson格式

86四、對流占主導地位時的數值方法引入人工擴散量方法8687四、對流占主導地位時的數值方法引入人工擴散量方法

為評價對流項引起的誤差，將上式各項歸并后得當人工擴散項與彌散權重

采用如下形式，對流項取非中心權而引起的誤差可消除：87四、對流占主導地位時的數值方法引入人工擴散量方法8788四、對流占主導地位時的數值方法引入人工擴散量方法

只有當時，方可導出對稱正定系數矩陣，上兩式可寫成

（6-116）（6-117）88四、對流占主導地位時的數值方法引入人工擴散量方法（6-18889四、對流占主導地位時的數值方法引入人工擴散量方法

（6-116）所定義的人工擴散是一種有效校正。

（6-117）彌散項取隱式對數值彌散沒有影響，卻提高求解精度。89四、對流占主導地位時的數值方法引入人工擴散量方法89地下水流數值模擬方法第六章

水動力彌散方程的數值解法中國地質大學環境學院2019春地下水流數值模擬方法第六章水動力彌散方程的數值解法中國地質9091一、有限差分法基本思想：

將空間劃分成若干網格，把時間分成若干小時段，每一個網格中心點（格點）處的未知變量（H或C）視為該網格平均值。利用差商代替微商。基本步驟：

（1）剖分（2）建立差分方程組（3）求解2一、有限差分法基本思想：9192一、有限差分法-導數的有限差分近似圖中，去x軸上任意一點i，其坐標為

在改點左右相聚為處分別取（i-1）和（i+1），其坐標分別為和以i為中心，泰勒展開C(x)3一、有限差分法-導數的有限差分近似圖中，去x軸上任意一點i9293一、有限差分法-導數的有限差分近似整理并略去余項（6-1）-（6-2），再除以略去余項（6-1）-（6-2），再除以略去余項：分別一階導數的向前差分、向后差分及中心差分二階中心差分4一、有限差分法-導數的有限差分近似整理并略去余項分別一階導9394一、有限差分法-一維水動力彌散的差分解法一維水動力彌散方程

（6-7）（1）顯格式式中僅有一個未知數，解得

式（6-7）中的對流項取中心差分5一、有限差分法-一維水動力彌散的差分解法一維水動力彌散方程9495

可以證明，穩定性準則要求

即

即。（1）（2）。由條件（2），格距要求很小；由（2）可知，鑒于較小，導致不能太大，將（2）代入（1）式中，得到

6

可以證明，穩定性準則要求

即9596若對流項改為向后差分

解得

穩定性要求不難看出，穩定性限制比對流項取中心差分有所改善。7若對流項改為向后差分解得9697（2）隱式格式

整理后得到

隱格式是無條件穩定的。8（2）隱式格式9798（2）Grank-Nicolson格式

整理后得到

取隱式和顯示的平均，即Grank-Nicolson格式9（2）Grank-Nicolson格式9899一、有限差分法-二維水動力彌散的差分解法二維水動力彌散方程

（4-56）（1）顯格式式中僅有一個未知數式（4-56）中的對流項取中心差分10一、有限差分法-二維水動力彌散的差分解法二維水動力彌散方99100一、有限差分法-二維水動力彌散的差分解法化簡后，有涉及以（i，j）為中心的5個網格點在tn時刻的已知濃度11一、有限差分法-二維水動力彌散的差分解法化簡后，有涉及以100101一、有限差分法-二維水動力彌散的差分解法（2）隱格式式（4-56）中右端的對流項取中心差分，右端個C的時階均取n+1水平12一、有限差分法-二維水動力彌散的差分解法（2）隱格式式（101102一、有限差分法-二維水動力彌散的差分解法（2）隱格式整理后收斂且無條件穩定涉及以（i，j）為中心的5個網格點在tn+1時刻的未知濃度13一、有限差分法-二維水動力彌散的差分解法（2）隱格式整理102103一、有限差分法-二維水動力彌散的差分解法（3）Grank-Nicolson格式將隱式格式的兩式相加除以214一、有限差分法-二維水動力彌散的差分解法（3）Grank103104一、有限差分法-二維水動力彌散的差分解法（3）Grank-Nicolson格式整理得15一、有限差分法-二維水動力彌散的差分解法（3）Grank104105一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱式法（ADI法）分兩次對三對角矩陣求逆，將一個Δt分成兩個Δt/2計算第一個半時間步，對x方向的偏導數采用隱式差分，對y方向的偏導數采用顯示差分。16一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱105106一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱式法（ADI法）整理得第二個半時間步，對y方向的偏導數采用隱式差分，對x方向的偏導數采用顯示差分。17一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱106107一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱式法（ADI法）整理得收斂且無條件穩定18一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱107108一、有限差分法-求解差分方程的計算機程序舉例算例

見教材P81-8519一、有限差分法-求解差分方程的計算機程序舉例算例108109二、有限單元法-伽遼金有限單元法伽遼金法

對均勻多孔介質，一維穩定流場中二維水動力彌散方程

取x軸方向與地下水流向相同，記伽遼金法是尋找一個級數形式的試探函數作微分方程的近似解，并使其滿足邊界條件

（6-26）20二、有限單元法-伽遼金有限單元法伽遼金法（6-26）109110二、有限單元法-伽遼金有限單元法上式一般不會滿足方程，因為僅是近似解，得到一個剩余誤差函數R(x，y)，在平均意義下迫使誤差為0

我們加以權，令剩余的加權積分為0，W是一組權函數。加權剩余法（6-29）21二、有限單元法-伽遼金有限單元法上式一般不會滿足方程，因110111二、有限單元法-伽遼金有限單元法在整個研究區D上，基函數NL(x,y)分段定義(1)函數區域連續（2）一階導數不連續（3）二階導不容易確定故采用分步積分的方法使二階導數降階22二、有限單元法-伽遼金有限單元法在整個研究區D上，基函數111112二、有限單元法-伽遼金有限單元法對一維積分整理后分步積分推廣到二維的情況（6-32）代入（6-29）23二、有限單元法-伽遼金有限單元法對一維積分（6-32）代112113二、有限單元法-伽遼金有限單元法有限單元剖分與基函數（1）三角形單元見《地下水流動問題數值模擬》（2）矩形單元將區域記為D，邊界記為B。要求各單元均質、等厚，即T、μ為常數。結點（內結點、邊界結點）

24二、有限單元法-伽遼金有限單元法有限單元剖分與基函數113114二、有限單元法-伽遼金有限單元法構造單元函數NL基函數取為“雙線性插值”

基函數NL（x，y）形狀如同一頂高等于1、有4條直線斜邊和4條下凹型曲邊的尖頂斗笠25二、有限單元法-伽遼金有限單元法構造單元函數NL基函數N114115二、有限單元法-伽遼金有限單元法子區DL以結點L為其共同結點的所有矩形單元（<=4）基函數NL僅在子區上不為0，在非DL上均為0,屬于子區DL單元e的單元基函數NL，用NeL表示（>0）26二、有限單元法-伽遼金有限單元法子區DL115116二、有限單元法-伽遼金有限單元法典型矩形單元該單元中心位于坐標原點處，且4條邊分別平行x，y軸,結點從左下角開始按逆時針方向編號i,j,m,n。沿x方向長2a，y方向寬2b。27二、有限單元法-伽遼金有限單元法典型矩形單元116117二、有限單元法-伽遼金有限單元法按上述要求所構造的NeL形式為對典型單元，令

則28二、有限單元法-伽遼金有限單元法按上述要求所構造的NeL117118二、有限單元法-伽遼金有限單元法其中即（6-34）29二、有限單元法-伽遼金有限單元法其中即（6-34）118119二、有限單元法-伽遼金有限單元法Nei(x,y)在結點i處為1，其它3處為0根據上式，有平行x或y方向上，Nei(x,y)線性變化故單元基函數為雙線性插值函數若結點坐標用表示結點濃度用表示30二、有限單元法-伽遼金有限單元法Nei(x,y)在結點i119120二、有限單元法-伽遼金有限單元法結點濃度CL(t)和單元基函數NeL(x,y)來確定單元內的近似解，即對于區域D，結點L的基函數在子區DL內各單元分塊確定。稱為矩陣單元e上的基函數（6-36）31二、有限單元法-伽遼金有限單元法結點濃度CL(t)和單元120121二、有限單元法-伽遼金有限單元法伽遼金有限單元法

左端積分范圍為D，由于在上NL=0,改寫在子區DL的積分，對矩形單元逐個積分、求和。設子區DL上有neL個單元，有（6-39）32二、有限單元法-伽遼金有限單元法伽遼金有限單元法（6-3121122二、有限單元法-伽遼金有限單元法由（6-36）得由（6-34）得再由（6-36）得（6-40）（6-42）33二、有限單元法-伽遼金有限單元法由（6-36）得由（6-122123二、有限單元法-伽遼金有限單元法將（6-40）（6-42）代入（6-39）得34二、有限單元法-伽遼金有限單元法將（6-40）（6-4123124二、有限單元法-伽遼金有限單元法（6-44）35二、有限單元法-伽遼金有限單元法（6-44）124125二、有限單元法-伽遼金有限單元法用矩陣形式表示一階微分方程組將一階導數用差分形式表示按矩陣符號有若﹛C﹜在新時間水平t+Δt上取值，則36二、有限單元法-伽遼金有限單元法用矩陣形式表示一階微分方125126二、有限單元法-伽遼金有限單元法（1）給定初始條件和邊界條件，可求（6-46）各個系數矩陣，于是可解得第一個末時刻各結點的濃度；（2）以此濃度為已知濃度，計算新的系數矩陣；（3）再次求解，得到第二個時間步長末時刻結點濃度（6-46）濃度可在之間任意位置近似表示37二、有限單元法-伽遼金有限單元法（1）給定初始條件和邊界126127二、有限單元法-伽遼金有限單元法按單元順序計算系數矩陣每個矩陣單元涉及4個結點i，j，m，n，涉及4個方程單元e的作用，可用4行與4列的矩陣來表示，稱為單元矩陣38二、有限單元法-伽遼金有限單元法按單元順序計算系數矩陣127128二、有限單元法-伽遼金有限單元法算例：單元彌散矩陣與整體彌散矩陣的關系設研究區劃分為8個單元，16個結點，四周為隔水邊界編號如下：39二、有限單元法-伽遼金有限單元法算例：單元彌散矩陣與整體128129二、有限單元法-伽遼金有限單元法算例：單元彌散矩陣與整體彌散矩陣的關系計算每個單元彌散矩陣，按雙下標已知，放置單元彌散矩陣并疊加，得到總體彌散矩陣。40二、有限單元法-伽遼金有限單元法算例：單元彌散矩陣與整體129130二、有限單元法-伽遼金有限單元法從（6-44）知將坐標原點移到矩形單元中心有41二、有限單元法-伽遼金有限單元法從（6-44）知將坐標原130131二、有限單元法-伽遼金有限單元法通過雙下標一致的元素求和，得到總體彌散矩陣諸元素42二、有限單元法-伽遼金有限單元法通過雙下標一致的元素求和131132二、有限單元法-伽遼金有限單元法積分求解可用高斯求積方法，詳見P97-10543二、有限單元法-伽遼金有限單元法積分求解可用高斯求積方法132133二、有限單元法-伽遼金有限單元法44二、有限單元法-伽遼金有限單元法133134二、有限單元法-伽遼金有限單元法算例：二維水動力彌散的伽遼金法計算機程序

數學模型寫成剖分成9個單元共20個結點，空間步長取Δx=Δy=10m，時間步長Δt=5d，滲透流速u=0.01m/d具體代碼見P106-11145二、有限單元法-伽遼金有限單元法算例：二維水動力彌散的伽134135三、過程與數值彌散過量現象

如一維流動水動力彌散中單位濃度梯度函數注入的情況。

對比解析解與數值解C-x曲線，在濃度封面附近數值計算的濃度超過最大濃度值1和小于最小濃度0，失去物理意義。數值彌散

若純對流方程采用有限差分逼近，其結果呈現似有彌散的存在。46三、過程與數值彌散過量現象數值彌散135136三、過程與數值彌散47三、過程與數值彌散136137三、過程與數值彌散過量現象-解析

由于時間步長與空間尺度匹配不佳，因而含水層在數值上不能“吸收”住污染物所引起。解決辦法：

將時間增量選擇為幾何級數的子項；

做試驗校正時間步長和差分網格；48三、過程與數值彌散過量現象-解析137138三、過程與數值彌散數值彌散-解析

由截斷誤差引起。

對流項一階空間導數采用向后差分逼近，濃度Ci-1泰勒展開變形得

向后差分的截斷誤差主部分為相當于彌散系數的彌散項49三、過程與數值彌散數值彌散-解析相當于彌散系數的彌散項138139三、過程與數值彌散數值彌散-解析

由故純對流方程有限差分近似表示后，結果如同對流-彌散方程的結果，記用數值彌散系數與物理彌散系數的比值來評價數值彌散對物理彌散的干涉程度此條件下的數值彌散系數50三、過程與數值彌散數值彌散-解析此條件下的數值彌散系數139140三、過程與數值彌散數值彌散-解析

定義為Peclet數，時，彌散占優；

時，對流占優。速度u與物理彌散系數DL比越大，越大。對于小數，數值解與精確解非常接近，對于大數，數值解過渡帶變寬，濃度鋒面變平緩，在鋒面的前后還出現解的振動—過量現象。控制Δx的大小從而減少數值彌散51三、過程與數值彌散數值彌散-解析控制Δx的大小從而減少數140141三、過程與數值彌散對純對流方程

52三、過程與數值彌散對純對流方程141142三、過程與數值彌散對純對流方程

結果也造成數值彌散此時，總的數值彌散是兩者之和53三、過程與數值彌散對純對流方程結果也造成數值彌散此時，總142143四、對流占主導地位時的數值方法上游加權法

對流項的有限差分化過程中乘上一個適當的權因子，則波動和過量得到改善或消除一維水動力彌散中對流項取中心隱式差分若相鄰格點i-1或（和）i+1的濃度越大，則格點i的濃度也越大，反映到差分方程上，系數Gi-1和Gi+1不得與Gi同號，上式要求

不能滿足時，對流占優勢54四、對流占主導地位時的數值方法上游加權法對流項的有限差分143144四、對流占主導地位時的數值方法上游加權法

假定DL和u為常數，取等格距差分網格。一維對流-彌散方程在積分，有對流-彌散通量用有限差分近似表示

55四、對流占主導地位時的數值方法上游加權法假定DL和u為常144145四、對流占主導地位時的數值方法上游加權法

ω為權因子；當ω=1/2時，對流項相當于中心差分；當ω=1時，對流項相當于取向后差分，

當ω=0時，對流項相當于向前差分。速度不同，

ω取值不同56四、對流占主導地位時的數值方法上游加權法145146四、對流占主導地位時的數值方法上游加權法

上游濃度權因子大于1/2。物理意義：加強上游格點濃度值在對流項中的作用可通過濃度的線性（或二次）插值具體確定在Δt期間C在格邊上隨時間的變化，從而確定ω值。57四、對流占主導地位時的數值方法上游加權法上游濃度權因子大146147四、對流占主導地位時的數值方法

仍用Gi-1、Gi和Gi+1分別表示上式左端第一、二、三項括號內系數。則采用上游加權法有限差分格式，可消除過量和波動，但截斷誤差增大，數值彌散反而增大。中心差分格式，截斷誤差為二階上游加權有限差分截斷誤差為一階58四、對流占主導地位時的數值方法仍用Gi-1、Gi147148四、對流占主導地位時的數值方法

上游權因子的物理實質：對流項濃度C在某斷面，例如斷面的中心差分意味著對流項的濃度在Δt時段內均以i與i-1兩格點的平均濃度，即通過斷面；實際為：Δt時段的初始時刻是該平均濃度值，隨后在對流的作用下，該斷面的濃度顯然會向上游格點的濃度Ci-1(u>0時)的變化。59四、對流占主導地位時的數值方法上游權因子的物理實148149四、對流占主導地位時的數值方法特征值方法

令源匯項為0，展開對流-彌散方程

60四、對流占主導地位時的數值方法特征值方法149150四、對流占主導地位時的數值方法特征值方法

假定流體密度值ρ不變（不可壓縮），含水層不可壓縮，則，（6-75）簡化為基本思路：物質導數表示一個沿軌跡（特征曲線）

運動的研究對象濃度變化率。由彌散作用引起。若存在運會想，還將受到沿軌跡源與匯、化學、生物化學反應、吸附與解析等。對流項用沿特征線的示蹤質點的運動描述

（6-76）61四、對流占主導地位時的數值方法特征值方法（6-76）150151四、對流占主導地位時的數值方法特征值方法

特征點均勻分布研究區，賦一個初始濃度。

含水層中未受污染部分特征點的初始濃度為零，特征點濃度眼軌跡的瞬時濃度變化，借助矩形網格差分格式計算。62四、對流占主導地位時的數值方法特征值方法151152四、對流占主導地位時的數值方法63四、對流占主導地位時的數值方法152153四、對流占主導地位時的數值方法模型精確度取決于網格距大小。可以在同一網格上計算速度場。如果采用較少的網格距，可以改善流場的描述。

坐標為的特征點所在的單元，可用下標i,j確定

精確度受到每一網格內所選擇的最初的特征點數目影響。64四、對流占主導地位時的數值方法模型精確度153154四、對流占主導地位時的數值方法缺點：（1）實際程序冗長，需記錄特征點的蹤跡；

（2）抽水點或出流邊界的特征點需要消除，入滲邊界須有特征點生成，不透水邊界特征點須反射回來；

（3）要通過消去特征點來避免滯留點附近特征點堆積；（4）特征點耗盡的單元，必須產生新的特征點。65四、對流占主導地位時的數值方法缺點：154155四、對流占主導地位時的數值方法動坐標系方法與網格變形方法能消除過量現象，還保持較陡的濃度鋒面。

對式

設坐標變換

得，

取空間步長，則

66四、對流占主導地位時的數值方法動坐標系方法與網格變形方法155156四、對流占主導地位時的數值方法動坐標系方法與網格變形方法

在固定坐標系中取格距為，有

動坐標系中格點運動與固定坐標系格點運動重合。

設固定坐標系中格點xj在tn+1時刻與動坐標系中Xi重合，有

動坐標系下濃度結點的表達式67四、對流占主導地位時的數值方法動坐標系方法與網格變形方法156157四、對流占主導地位時的數值方法動坐標系方法與網格變形方法

式變成

對半無限長砂柱，當DL=0時68四、對流占主導地位時的數值方法動坐標系方法與網格變形方法157158四、對流占主導地位時的數值方法網格變形方法

設tn時刻有限元結點在空間中的分布是對應濃度分布為

69四、對流占主導地位時的數值方法網格變形方法158159四、對流占主導地位時的數值方法網格變形方法

主要步驟：結點位置隨時間變化，基函數及有限元方程的系數都依賴于時間，每推進一個Δt都需要重新計算。二維問題還須進行畸形判斷

70四、對流占主導地位時的數值方法網格變形方法159160四、對流占主導地位時的數值方法結合動坐標的網格變形方法

一維對流-彌散在某個動點在t時刻位置x可寫成

其中，x0是動點的初始位置，有取時間的一階偏導數對流-彌散方程寫成

若控制動點速度dx/dt并讓它接近流速u，則可轉換成彌散為主的的方程式71四、對流占主導地位時的數值方法結合動坐標的網格變形方法若160161四、對流占主導地位時的數值方法若使用伽遼金法，基函數因依賴于結點位置而成為時間的函數。若用表示隨時間變化的基函數，則對任意t，有

故72四、對流占主導地位時的數值方法若使用伽遼金法，基函數因依161162四、對流占主