溶質(zhì)運(yùn)移理論水動(dòng)力彌散方程的值解法-課件

地下水流數(shù)值模擬方法第六章

水動(dòng)力彌散方程的數(shù)值解法中國(guó)地質(zhì)大學(xué)環(huán)境學(xué)院2019春地下水流數(shù)值模擬方法第六章水動(dòng)力彌散方程的數(shù)值解法中國(guó)地質(zhì)12一、有限差分法基本思想：

將空間劃分成若干網(wǎng)格，把時(shí)間分成若干小時(shí)段，每一個(gè)網(wǎng)格中心點(diǎn)（格點(diǎn)）處的未知變量（H或C）視為該網(wǎng)格平均值。利用差商代替微商?；静襟E：

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

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

（6-7）（1）顯格式式中僅有一個(gè)未知數(shù)，解得

式（6-7）中的對(duì)流項(xiàng)取中心差分5一、有限差分法-一維水動(dòng)力彌散的差分解法一維水動(dòng)力彌散方程56

可以證明，穩(wěn)定性準(zhǔn)則要求

即

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

6

可以證明，穩(wěn)定性準(zhǔn)則要求

即67若對(duì)流項(xiàng)改為向后差分

解得

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

整理后得到

隱格式是無(wú)條件穩(wěn)定的。8（2）隱式格式89（2）Grank-Nicolson格式

整理后得到

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

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

見(jiàn)教材P81-8519一、有限差分法-求解差分方程的計(jì)算機(jī)程序舉例算例1920二、有限單元法-伽遼金有限單元法伽遼金法

對(duì)均勻多孔介質(zhì)，一維穩(wěn)定流場(chǎng)中二維水動(dòng)力彌散方程

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

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

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

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

基函數(shù)NL（x，y）形狀如同一頂高等于1、有4條直線斜邊和4條下凹型曲邊的尖頂斗笠25二、有限單元法-伽遼金有限單元法構(gòu)造單元函數(shù)NL基函數(shù)N2526二、有限單元法-伽遼金有限單元法子區(qū)DL以結(jié)點(diǎn)L為其共同結(jié)點(diǎn)的所有矩形單元（<=4）基函數(shù)NL僅在子區(qū)上不為0，在非DL上均為0,屬于子區(qū)DL單元e的單元基函數(shù)NL，用NeL表示（>0）26二、有限單元法-伽遼金有限單元法子區(qū)DL2627二、有限單元法-伽遼金有限單元法典型矩形單元該單元中心位于坐標(biāo)原點(diǎn)處，且4條邊分別平行x，y軸,結(jié)點(diǎn)從左下角開(kāi)始按逆時(shí)針?lè)较蚓幪?hào)i,j,m,n。沿x方向長(zhǎng)2a，y方向?qū)?b。27二、有限單元法-伽遼金有限單元法典型矩形單元2728二、有限單元法-伽遼金有限單元法按上述要求所構(gòu)造的NeL形式為對(duì)典型單元，令

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

左端積分范圍為D，由于在上NL=0,改寫在子區(qū)DL的積分，對(duì)矩形單元逐個(gè)積分、求和。設(shè)子區(qū)DL上有neL個(gè)單元，有（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二、有限單元法-伽遼金有限單元法用矩陣形式表示一階微分方程組將一階導(dǎo)數(shù)用差分形式表示按矩陣符號(hào)有若﹛C﹜在新時(shí)間水平t+Δt上取值，則36二、有限單元法-伽遼金有限單元法用矩陣形式表示一階微分方3637二、有限單元法-伽遼金有限單元法（1）給定初始條件和邊界條件，可求（6-46）各個(gè)系數(shù)矩陣，于是可解得第一個(gè)末時(shí)刻各結(jié)點(diǎn)的濃度；（2）以此濃度為已知濃度，計(jì)算新的系數(shù)矩陣；（3）再次求解，得到第二個(gè)時(shí)間步長(zhǎng)末時(shí)刻結(jié)點(diǎn)濃度（6-46）濃度可在之間任意位置近似表示37二、有限單元法-伽遼金有限單元法（1）給定初始條件和邊界3738二、有限單元法-伽遼金有限單元法按單元順序計(jì)算系數(shù)矩陣每個(gè)矩陣單元涉及4個(gè)結(jié)點(diǎn)i，j，m，n，涉及4個(gè)方程單元e的作用，可用4行與4列的矩陣來(lái)表示，稱為單元矩陣38二、有限單元法-伽遼金有限單元法按單元順序計(jì)算系數(shù)矩陣3839二、有限單元法-伽遼金有限單元法算例：?jiǎn)卧獜浬⒕仃嚺c整體彌散矩陣的關(guān)系設(shè)研究區(qū)劃分為8個(gè)單元，16個(gè)結(jié)點(diǎn)，四周為隔水邊界編號(hào)如下：39二、有限單元法-伽遼金有限單元法算例：?jiǎn)卧獜浬⒕仃嚺c整體3940二、有限單元法-伽遼金有限單元法算例：?jiǎn)卧獜浬⒕仃嚺c整體彌散矩陣的關(guān)系計(jì)算每個(gè)單元彌散矩陣，按雙下標(biāo)已知，放置單元彌散矩陣并疊加，得到總體彌散矩陣。40二、有限單元法-伽遼金有限單元法算例：?jiǎn)卧獜浬⒕仃嚺c整體4041二、有限單元法-伽遼金有限單元法從（6-44）知將坐標(biāo)原點(diǎn)移到矩形單元中心有41二、有限單元法-伽遼金有限單元法從（6-44）知將坐標(biāo)原4142二、有限單元法-伽遼金有限單元法通過(guò)雙下標(biāo)一致的元素求和，得到總體彌散矩陣諸元素42二、有限單元法-伽遼金有限單元法通過(guò)雙下標(biāo)一致的元素求和4243二、有限單元法-伽遼金有限單元法積分求解可用高斯求積方法，詳見(jiàn)P97-10543二、有限單元法-伽遼金有限單元法積分求解可用高斯求積方法4344二、有限單元法-伽遼金有限單元法44二、有限單元法-伽遼金有限單元法4445二、有限單元法-伽遼金有限單元法算例：二維水動(dòng)力彌散的伽遼金法計(jì)算機(jī)程序

數(shù)學(xué)模型寫成剖分成9個(gè)單元共20個(gè)結(jié)點(diǎn)，空間步長(zhǎng)取Δx=Δy=10m，時(shí)間步長(zhǎng)Δt=5d，滲透流速u=0.01m/d具體代碼見(jiàn)P106-11145二、有限單元法-伽遼金有限單元法算例：二維水動(dòng)力彌散的伽4546三、過(guò)程與數(shù)值彌散過(guò)量現(xiàn)象

如一維流動(dòng)水動(dòng)力彌散中單位濃度梯度函數(shù)注入的情況。

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

若純對(duì)流方程采用有限差分逼近，其結(jié)果呈現(xiàn)似有彌散的存在。46三、過(guò)程與數(shù)值彌散過(guò)量現(xiàn)象數(shù)值彌散4647三、過(guò)程與數(shù)值彌散47三、過(guò)程與數(shù)值彌散4748三、過(guò)程與數(shù)值彌散過(guò)量現(xiàn)象-解析

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

將時(shí)間增量選擇為幾何級(jí)數(shù)的子項(xiàng)；

做試驗(yàn)校正時(shí)間步長(zhǎng)和差分網(wǎng)格；48三、過(guò)程與數(shù)值彌散過(guò)量現(xiàn)象-解析4849三、過(guò)程與數(shù)值彌散數(shù)值彌散-解析

由截?cái)嗾`差引起。

對(duì)流項(xiàng)一階空間導(dǎo)數(shù)采用向后差分逼近，濃度Ci-1泰勒展開(kāi)變形得

向后差分的截?cái)嗾`差主部分為相當(dāng)于彌散系數(shù)的彌散項(xiàng)49三、過(guò)程與數(shù)值彌散數(shù)值彌散-解析相當(dāng)于彌散系數(shù)的彌散項(xiàng)4950三、過(guò)程與數(shù)值彌散數(shù)值彌散-解析

由故純對(duì)流方程有限差分近似表示后，結(jié)果如同對(duì)流-彌散方程的結(jié)果，記用數(shù)值彌散系數(shù)與物理彌散系數(shù)的比值來(lái)評(píng)價(jià)數(shù)值彌散對(duì)物理彌散的干涉程度此條件下的數(shù)值彌散系數(shù)50三、過(guò)程與數(shù)值彌散數(shù)值彌散-解析此條件下的數(shù)值彌散系數(shù)5051三、過(guò)程與數(shù)值彌散數(shù)值彌散-解析

定義為Peclet數(shù)，時(shí)，彌散占優(yōu)；

時(shí)，對(duì)流占優(yōu)。速度u與物理彌散系數(shù)DL比越大，越大。對(duì)于小數(shù)，數(shù)值解與精確解非常接近，對(duì)于大數(shù)，數(shù)值解過(guò)渡帶變寬，濃度鋒面變平緩，在鋒面的前后還出現(xiàn)解的振動(dòng)—過(guò)量現(xiàn)象?？刂痞的大小從而減少數(shù)值彌散51三、過(guò)程與數(shù)值彌散數(shù)值彌散-解析控制Δx的大小從而減少數(shù)5152三、過(guò)程與數(shù)值彌散對(duì)純對(duì)流方程

52三、過(guò)程與數(shù)值彌散對(duì)純對(duì)流方程5253三、過(guò)程與數(shù)值彌散對(duì)純對(duì)流方程

結(jié)果也造成數(shù)值彌散此時(shí)，總的數(shù)值彌散是兩者之和53三、過(guò)程與數(shù)值彌散對(duì)純對(duì)流方程結(jié)果也造成數(shù)值彌散此時(shí)，總5354四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法

對(duì)流項(xiàng)的有限差分化過(guò)程中乘上一個(gè)適當(dāng)?shù)臋?quán)因子，則波動(dòng)和過(guò)量得到改善或消除一維水動(dòng)力彌散中對(duì)流項(xiàng)取中心隱式差分若相鄰格點(diǎn)i-1或（和）i+1的濃度越大，則格點(diǎn)i的濃度也越大，反映到差分方程上，系數(shù)Gi-1和Gi+1不得與Gi同號(hào)，上式要求

不能滿足時(shí)，對(duì)流占優(yōu)勢(shì)54四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法對(duì)流項(xiàng)的有限差分5455四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法

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

55四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法假定DL和u為常5556四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法

ω為權(quán)因子；當(dāng)ω=1/2時(shí)，對(duì)流項(xiàng)相當(dāng)于中心差分；當(dāng)ω=1時(shí)，對(duì)流項(xiàng)相當(dāng)于取向后差分，

當(dāng)ω=0時(shí)，對(duì)流項(xiàng)相當(dāng)于向前差分。速度不同，

ω取值不同56四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法5657四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法

上游濃度權(quán)因子大于1/2。物理意義：加強(qiáng)上游格點(diǎn)濃度值在對(duì)流項(xiàng)中的作用可通過(guò)濃度的線性（或二次）插值具體確定在Δt期間C在格邊上隨時(shí)間的變化，從而確定ω值。57四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法上游濃度權(quán)因子大5758四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法

仍用Gi-1、Gi和Gi+1分別表示上式左端第一、二、三項(xiàng)括號(hào)內(nèi)系數(shù)。則采用上游加權(quán)法有限差分格式，可消除過(guò)量和波動(dòng)，但截?cái)嗾`差增大，數(shù)值彌散反而增大。中心差分格式，截?cái)嗾`差為二階上游加權(quán)有限差分截?cái)嗾`差為一階58四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法仍用Gi-1、Gi5859四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法

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

令源匯項(xiàng)為0，展開(kāi)對(duì)流-彌散方程

60四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法特征值方法6061四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法特征值方法

假定流體密度值ρ不變（不可壓縮），含水層不可壓縮，則，（6-75）簡(jiǎn)化為基本思路：物質(zhì)導(dǎo)數(shù)表示一個(gè)沿軌跡（特征曲線）

運(yùn)動(dòng)的研究對(duì)象濃度變化率。由彌散作用引起。若存在運(yùn)會(huì)想，還將受到沿軌跡源與匯、化學(xué)、生物化學(xué)反應(yīng)、吸附與解析等。對(duì)流項(xiàng)用沿特征線的示蹤質(zhì)點(diǎn)的運(yùn)動(dòng)描述

（6-76）61四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法特征值方法（6-76）6162四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法特征值方法

特征點(diǎn)均勻分布研究區(qū)，賦一個(gè)初始濃度。

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

坐標(biāo)為的特征點(diǎn)所在的單元，可用下標(biāo)i,j確定

精確度受到每一網(wǎng)格內(nèi)所選擇的最初的特征點(diǎn)數(shù)目影響。64四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法模型精確度6465四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法缺點(diǎn)：（1）實(shí)際程序冗長(zhǎng)，需記錄特征點(diǎn)的蹤跡；

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

（3）要通過(guò)消去特征點(diǎn)來(lái)避免滯留點(diǎn)附近特征點(diǎn)堆積；（4）特征點(diǎn)耗盡的單元，必須產(chǎn)生新的特征點(diǎn)。65四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法缺點(diǎn)：6566四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法動(dòng)坐標(biāo)系方法與網(wǎng)格變形方法能消除過(guò)量現(xiàn)象，還保持較陡的濃度鋒面。

對(duì)式

設(shè)坐標(biāo)變換

得，

取空間步長(zhǎng)，則

66四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法動(dòng)坐標(biāo)系方法與網(wǎng)格變形方法6667四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法動(dòng)坐標(biāo)系方法與網(wǎng)格變形方法

在固定坐標(biāo)系中取格距為，有

動(dòng)坐標(biāo)系中格點(diǎn)運(yùn)動(dòng)與固定坐標(biāo)系格點(diǎn)運(yùn)動(dòng)重合。

設(shè)固定坐標(biāo)系中格點(diǎn)xj在tn+1時(shí)刻與動(dòng)坐標(biāo)系中Xi重合，有

動(dòng)坐標(biāo)系下濃度結(jié)點(diǎn)的表達(dá)式67四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法動(dòng)坐標(biāo)系方法與網(wǎng)格變形方法6768四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法動(dòng)坐標(biāo)系方法與網(wǎng)格變形方法

式變成

對(duì)半無(wú)限長(zhǎng)砂柱，當(dāng)DL=0時(shí)68四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法動(dòng)坐標(biāo)系方法與網(wǎng)格變形方法6869四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法網(wǎng)格變形方法

設(shè)tn時(shí)刻有限元結(jié)點(diǎn)在空間中的分布是對(duì)應(yīng)濃度分布為

69四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法網(wǎng)格變形方法6970四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法網(wǎng)格變形方法

主要步驟：結(jié)點(diǎn)位置隨時(shí)間變化，基函數(shù)及有限元方程的系數(shù)都依賴于時(shí)間，每推進(jìn)一個(gè)Δt都需要重新計(jì)算。二維問(wèn)題還須進(jìn)行畸形判斷

70四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法網(wǎng)格變形方法7071四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法結(jié)合動(dòng)坐標(biāo)的網(wǎng)格變形方法

一維對(duì)流-彌散在某個(gè)動(dòng)點(diǎn)在t時(shí)刻位置x可寫成

其中，x0是動(dòng)點(diǎn)的初始位置，有取時(shí)間的一階偏導(dǎo)數(shù)對(duì)流-彌散方程寫成

若控制動(dòng)點(diǎn)速度dx/dt并讓它接近流速u，則可轉(zhuǎn)換成彌散為主的的方程式71四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法結(jié)合動(dòng)坐標(biāo)的網(wǎng)格變形方法若7172四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法若使用伽遼金法，基函數(shù)因依賴于結(jié)點(diǎn)位置而成為時(shí)間的函數(shù)。若用表示隨時(shí)間變化的基函數(shù)，則對(duì)任意t，有

故72四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法若使用伽遼金法，基函數(shù)因依7273四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法

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

73四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法7374四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法-以一維運(yùn)移為例在x=0處，一個(gè)瞬時(shí)注入的、質(zhì)量為Δ

M的理想失蹤及，在t時(shí)刻濃度分布對(duì)一個(gè)固定時(shí)間t，看做正態(tài)分布74四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法-以一維運(yùn)移為例7475四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法-以一維運(yùn)移為例在時(shí)間t=0和位置x=0處，放置具有污染物質(zhì)量為ΔM/N的N個(gè)特征點(diǎn)，每一個(gè)特征點(diǎn)移動(dòng)距離x，而在時(shí)刻t到達(dá)它們各自位置，即Z是一個(gè)平均值為0，標(biāo)準(zhǔn)差為1的正態(tài)分布隨機(jī)變量

所得的頻率分布，通過(guò)一個(gè)標(biāo)準(zhǔn)化因子，有對(duì)有限數(shù)量的顆粒，可近似去頂C(x,t)，將空間變量劃分為長(zhǎng)度Δx的若干區(qū)間，區(qū)間中點(diǎn)濃度75四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法-以一維運(yùn)移為例7576四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法-以一維運(yùn)移為例若孔隙平均流速是時(shí)間和位置的函數(shù)，那么時(shí)間為0開(kāi)始的顆粒P在時(shí)間間隔Δt內(nèi)的隨機(jī)軌跡xp(t)為

76四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法-以一維運(yùn)移為例7677四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法-二維流場(chǎng)

首先假定流動(dòng)方向平行于x軸，得到Z和Z’是正態(tài)分布的隨機(jī)變量的兩個(gè)值。對(duì)任意方向的流動(dòng)，需考慮彌散張量特性。77四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法-二維流場(chǎng)7778四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法-二維流場(chǎng)

首先假定流動(dòng)方向平行于x軸，得到Z和Z’是正態(tài)分布的隨機(jī)變量的兩個(gè)值。對(duì)任意方向的流動(dòng)，需考慮彌散張量特性。對(duì)流運(yùn)動(dòng)可確定在水流流動(dòng)方向及其正交方向的彌散運(yùn)動(dòng)78四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法-二維流場(chǎng)對(duì)流運(yùn)7879四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法

79四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法7980四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法

數(shù)值結(jié)果的質(zhì)量主要取決于特征點(diǎn)數(shù)目N，要考慮網(wǎng)格的選擇。大網(wǎng)格間距會(huì)產(chǎn)生極平均的結(jié)果，小網(wǎng)格間距會(huì)產(chǎn)生十分粗糙的分布。

在穩(wěn)定流動(dòng)和一個(gè)源或數(shù)個(gè)源都有固定強(qiáng)度的情況下，可以節(jié)省很多特征點(diǎn)。80四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法8081818182四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法

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

不會(huì)出現(xiàn)數(shù)值彌散，但在污染帶邊緣的得到滿意的計(jì)算結(jié)果需要大量的特征點(diǎn)。

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

靈敏度較小的參數(shù)的變化引起的濃度分布變化會(huì)被隨機(jī)性濃度的變化覆蓋。

83四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法隨機(jī)步行法-特點(diǎn)8384四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法引入人工擴(kuò)散量方法

對(duì)流-彌散方程簡(jiǎn)寫成簡(jiǎn)寫成腳標(biāo)注釋的形式引入人工擴(kuò)散系數(shù)并賦權(quán)得

（6-107）84四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法引入人工擴(kuò)散量方法（6-18485四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法引入人工擴(kuò)散量方法

當(dāng)θu=0時(shí)，導(dǎo)出對(duì)稱的正定矩陣，若不校正，格式是不穩(wěn)定的。在時(shí)間步長(zhǎng)中點(diǎn)泰勒展開(kāi)

（6-109）85四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法引入人工擴(kuò)散量方法（6-18586四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法引入人工擴(kuò)散量方法

（6-107）對(duì)時(shí)間求偏導(dǎo)并代入（6-109）得在時(shí)，是具有二階截?cái)嗾`差的Grank-Nicolson格式

86四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法引入人工擴(kuò)散量方法8687四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法引入人工擴(kuò)散量方法

為評(píng)價(jià)對(duì)流項(xiàng)引起的誤差，將上式各項(xiàng)歸并后得當(dāng)人工擴(kuò)散項(xiàng)與彌散權(quán)重

采用如下形式，對(duì)流項(xiàng)取非中心權(quán)而引起的誤差可消除：87四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法引入人工擴(kuò)散量方法8788四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法引入人工擴(kuò)散量方法

只有當(dāng)時(shí)，方可導(dǎo)出對(duì)稱正定系數(shù)矩陣，上兩式可寫成

（6-116）（6-117）88四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法引入人工擴(kuò)散量方法（6-18889四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法引入人工擴(kuò)散量方法

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

（6-117）彌散項(xiàng)取隱式對(duì)數(shù)值彌散沒(méi)有影響，卻提高求解精度。89四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法引入人工擴(kuò)散量方法89地下水流數(shù)值模擬方法第六章

水動(dòng)力彌散方程的數(shù)值解法中國(guó)地質(zhì)大學(xué)環(huán)境學(xué)院2019春地下水流數(shù)值模擬方法第六章水動(dòng)力彌散方程的數(shù)值解法中國(guó)地質(zhì)9091一、有限差分法基本思想：

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

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

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

（6-7）（1）顯格式式中僅有一個(gè)未知數(shù)，解得

式（6-7）中的對(duì)流項(xiàng)取中心差分5一、有限差分法-一維水動(dòng)力彌散的差分解法一維水動(dòng)力彌散方程9495

可以證明，穩(wěn)定性準(zhǔn)則要求

即

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

6

可以證明，穩(wěn)定性準(zhǔn)則要求

即9596若對(duì)流項(xiàng)改為向后差分

解得

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

整理后得到

隱格式是無(wú)條件穩(wěn)定的。8（2）隱式格式9798（2）Grank-Nicolson格式

整理后得到

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

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

見(jiàn)教材P81-8519一、有限差分法-求解差分方程的計(jì)算機(jī)程序舉例算例108109二、有限單元法-伽遼金有限單元法伽遼金法

對(duì)均勻多孔介質(zhì)，一維穩(wěn)定流場(chǎng)中二維水動(dòng)力彌散方程

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

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

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

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

基函數(shù)NL（x，y）形狀如同一頂高等于1、有4條直線斜邊和4條下凹型曲邊的尖頂斗笠25二、有限單元法-伽遼金有限單元法構(gòu)造單元函數(shù)NL基函數(shù)N114115二、有限單元法-伽遼金有限單元法子區(qū)DL以結(jié)點(diǎn)L為其共同結(jié)點(diǎn)的所有矩形單元（<=4）基函數(shù)NL僅在子區(qū)上不為0，在非DL上均為0,屬于子區(qū)DL單元e的單元基函數(shù)NL，用NeL表示（>0）26二、有限單元法-伽遼金有限單元法子區(qū)DL115116二、有限單元法-伽遼金有限單元法典型矩形單元該單元中心位于坐標(biāo)原點(diǎn)處，且4條邊分別平行x，y軸,結(jié)點(diǎn)從左下角開(kāi)始按逆時(shí)針?lè)较蚓幪?hào)i,j,m,n。沿x方向長(zhǎng)2a，y方向?qū)?b。27二、有限單元法-伽遼金有限單元法典型矩形單元116117二、有限單元法-伽遼金有限單元法按上述要求所構(gòu)造的NeL形式為對(duì)典型單元，令

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

左端積分范圍為D，由于在上NL=0,改寫在子區(qū)DL的積分，對(duì)矩形單元逐個(gè)積分、求和。設(shè)子區(qū)DL上有neL個(gè)單元，有（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二、有限單元法-伽遼金有限單元法用矩陣形式表示一階微分方程組將一階導(dǎo)數(shù)用差分形式表示按矩陣符號(hào)有若﹛C﹜在新時(shí)間水平t+Δt上取值，則36二、有限單元法-伽遼金有限單元法用矩陣形式表示一階微分方125126二、有限單元法-伽遼金有限單元法（1）給定初始條件和邊界條件，可求（6-46）各個(gè)系數(shù)矩陣，于是可解得第一個(gè)末時(shí)刻各結(jié)點(diǎn)的濃度；（2）以此濃度為已知濃度，計(jì)算新的系數(shù)矩陣；（3）再次求解，得到第二個(gè)時(shí)間步長(zhǎng)末時(shí)刻結(jié)點(diǎn)濃度（6-46）濃度可在之間任意位置近似表示37二、有限單元法-伽遼金有限單元法（1）給定初始條件和邊界126127二、有限單元法-伽遼金有限單元法按單元順序計(jì)算系數(shù)矩陣每個(gè)矩陣單元涉及4個(gè)結(jié)點(diǎn)i，j，m，n，涉及4個(gè)方程單元e的作用，可用4行與4列的矩陣來(lái)表示，稱為單元矩陣38二、有限單元法-伽遼金有限單元法按單元順序計(jì)算系數(shù)矩陣127128二、有限單元法-伽遼金有限單元法算例：?jiǎn)卧獜浬⒕仃嚺c整體彌散矩陣的關(guān)系設(shè)研究區(qū)劃分為8個(gè)單元，16個(gè)結(jié)點(diǎn)，四周為隔水邊界編號(hào)如下：39二、有限單元法-伽遼金有限單元法算例：?jiǎn)卧獜浬⒕仃嚺c整體128129二、有限單元法-伽遼金有限單元法算例：?jiǎn)卧獜浬⒕仃嚺c整體彌散矩陣的關(guān)系計(jì)算每個(gè)單元彌散矩陣，按雙下標(biāo)已知，放置單元彌散矩陣并疊加，得到總體彌散矩陣。40二、有限單元法-伽遼金有限單元法算例：?jiǎn)卧獜浬⒕仃嚺c整體129130二、有限單元法-伽遼金有限單元法從（6-44）知將坐標(biāo)原點(diǎn)移到矩形單元中心有41二、有限單元法-伽遼金有限單元法從（6-44）知將坐標(biāo)原130131二、有限單元法-伽遼金有限單元法通過(guò)雙下標(biāo)一致的元素求和，得到總體彌散矩陣諸元素42二、有限單元法-伽遼金有限單元法通過(guò)雙下標(biāo)一致的元素求和131132二、有限單元法-伽遼金有限單元法積分求解可用高斯求積方法，詳見(jiàn)P97-10543二、有限單元法-伽遼金有限單元法積分求解可用高斯求積方法132133二、有限單元法-伽遼金有限單元法44二、有限單元法-伽遼金有限單元法133134二、有限單元法-伽遼金有限單元法算例：二維水動(dòng)力彌散的伽遼金法計(jì)算機(jī)程序

數(shù)學(xué)模型寫成剖分成9個(gè)單元共20個(gè)結(jié)點(diǎn)，空間步長(zhǎng)取Δx=Δy=10m，時(shí)間步長(zhǎng)Δt=5d，滲透流速u=0.01m/d具體代碼見(jiàn)P106-11145二、有限單元法-伽遼金有限單元法算例：二維水動(dòng)力彌散的伽134135三、過(guò)程與數(shù)值彌散過(guò)量現(xiàn)象

如一維流動(dòng)水動(dòng)力彌散中單位濃度梯度函數(shù)注入的情況。

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

若純對(duì)流方程采用有限差分逼近，其結(jié)果呈現(xiàn)似有彌散的存在。46三、過(guò)程與數(shù)值彌散過(guò)量現(xiàn)象數(shù)值彌散135136三、過(guò)程與數(shù)值彌散47三、過(guò)程與數(shù)值彌散136137三、過(guò)程與數(shù)值彌散過(guò)量現(xiàn)象-解析

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

將時(shí)間增量選擇為幾何級(jí)數(shù)的子項(xiàng)；

做試驗(yàn)校正時(shí)間步長(zhǎng)和差分網(wǎng)格；48三、過(guò)程與數(shù)值彌散過(guò)量現(xiàn)象-解析137138三、過(guò)程與數(shù)值彌散數(shù)值彌散-解析

由截?cái)嗾`差引起。

對(duì)流項(xiàng)一階空間導(dǎo)數(shù)采用向后差分逼近，濃度Ci-1泰勒展開(kāi)變形得

向后差分的截?cái)嗾`差主部分為相當(dāng)于彌散系數(shù)的彌散項(xiàng)49三、過(guò)程與數(shù)值彌散數(shù)值彌散-解析相當(dāng)于彌散系數(shù)的彌散項(xiàng)138139三、過(guò)程與數(shù)值彌散數(shù)值彌散-解析

由故純對(duì)流方程有限差分近似表示后，結(jié)果如同對(duì)流-彌散方程的結(jié)果，記用數(shù)值彌散系數(shù)與物理彌散系數(shù)的比值來(lái)評(píng)價(jià)數(shù)值彌散對(duì)物理彌散的干涉程度此條件下的數(shù)值彌散系數(shù)50三、過(guò)程與數(shù)值彌散數(shù)值彌散-解析此條件下的數(shù)值彌散系數(shù)139140三、過(guò)程與數(shù)值彌散數(shù)值彌散-解析

定義為Peclet數(shù)，時(shí)，彌散占優(yōu)；

時(shí)，對(duì)流占優(yōu)。速度u與物理彌散系數(shù)DL比越大，越大。對(duì)于小數(shù)，數(shù)值解與精確解非常接近，對(duì)于大數(shù)，數(shù)值解過(guò)渡帶變寬，濃度鋒面變平緩，在鋒面的前后還出現(xiàn)解的振動(dòng)—過(guò)量現(xiàn)象?？刂痞的大小從而減少數(shù)值彌散51三、過(guò)程與數(shù)值彌散數(shù)值彌散-解析控制Δx的大小從而減少數(shù)140141三、過(guò)程與數(shù)值彌散對(duì)純對(duì)流方程

52三、過(guò)程與數(shù)值彌散對(duì)純對(duì)流方程141142三、過(guò)程與數(shù)值彌散對(duì)純對(duì)流方程

結(jié)果也造成數(shù)值彌散此時(shí)，總的數(shù)值彌散是兩者之和53三、過(guò)程與數(shù)值彌散對(duì)純對(duì)流方程結(jié)果也造成數(shù)值彌散此時(shí)，總142143四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法

對(duì)流項(xiàng)的有限差分化過(guò)程中乘上一個(gè)適當(dāng)?shù)臋?quán)因子，則波動(dòng)和過(guò)量得到改善或消除一維水動(dòng)力彌散中對(duì)流項(xiàng)取中心隱式差分若相鄰格點(diǎn)i-1或（和）i+1的濃度越大，則格點(diǎn)i的濃度也越大，反映到差分方程上，系數(shù)Gi-1和Gi+1不得與Gi同號(hào)，上式要求

不能滿足時(shí)，對(duì)流占優(yōu)勢(shì)54四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法對(duì)流項(xiàng)的有限差分143144四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法

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

55四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法假定DL和u為常144145四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法

ω為權(quán)因子；當(dāng)ω=1/2時(shí)，對(duì)流項(xiàng)相當(dāng)于中心差分；當(dāng)ω=1時(shí)，對(duì)流項(xiàng)相當(dāng)于取向后差分，

當(dāng)ω=0時(shí)，對(duì)流項(xiàng)相當(dāng)于向前差分。速度不同，

ω取值不同56四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法145146四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法

上游濃度權(quán)因子大于1/2。物理意義：加強(qiáng)上游格點(diǎn)濃度值在對(duì)流項(xiàng)中的作用可通過(guò)濃度的線性（或二次）插值具體確定在Δt期間C在格邊上隨時(shí)間的變化，從而確定ω值。57四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法上游加權(quán)法上游濃度權(quán)因子大146147四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法

仍用Gi-1、Gi和Gi+1分別表示上式左端第一、二、三項(xiàng)括號(hào)內(nèi)系數(shù)。則采用上游加權(quán)法有限差分格式，可消除過(guò)量和波動(dòng)，但截?cái)嗾`差增大，數(shù)值彌散反而增大。中心差分格式，截?cái)嗾`差為二階上游加權(quán)有限差分截?cái)嗾`差為一階58四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法仍用Gi-1、Gi147148四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法

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

令源匯項(xiàng)為0，展開(kāi)對(duì)流-彌散方程

60四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法特征值方法149150四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法特征值方法

假定流體密度值ρ不變（不可壓縮），含水層不可壓縮，則，（6-75）簡(jiǎn)化為基本思路：物質(zhì)導(dǎo)數(shù)表示一個(gè)沿軌跡（特征曲線）

運(yùn)動(dòng)的研究對(duì)象濃度變化率。由彌散作用引起。若存在運(yùn)會(huì)想，還將受到沿軌跡源與匯、化學(xué)、生物化學(xué)反應(yīng)、吸附與解析等。對(duì)流項(xiàng)用沿特征線的示蹤質(zhì)點(diǎn)的運(yùn)動(dòng)描述

（6-76）61四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法特征值方法（6-76）150151四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法特征值方法

特征點(diǎn)均勻分布研究區(qū)，賦一個(gè)初始濃度。

含水層中未受污染部分特征點(diǎn)的初始濃度為零，特征點(diǎn)濃度眼軌跡的瞬時(shí)濃度變化，借助矩形網(wǎng)格差分格式計(jì)算。62四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法特征值方法151152四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法63四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法152153四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法模型精確度取決于網(wǎng)格距大小?？梢栽谕痪W(wǎng)格上計(jì)算速度場(chǎng)。如果采用較少的網(wǎng)格距，可以改善流場(chǎng)的描述。

坐標(biāo)為的特征點(diǎn)所在的單元，可用下標(biāo)i,j確定

精確度受到每一網(wǎng)格內(nèi)所選擇的最初的特征點(diǎn)數(shù)目影響。64四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法模型精確度153154四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法缺點(diǎn)：（1）實(shí)際程序冗長(zhǎng)，需記錄特征點(diǎn)的蹤跡；

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

（3）要通過(guò)消去特征點(diǎn)來(lái)避免滯留點(diǎn)附近特征點(diǎn)堆積；（4）特征點(diǎn)耗盡的單元，必須產(chǎn)生新的特征點(diǎn)。65四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法缺點(diǎn)：154155四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法動(dòng)坐標(biāo)系方法與網(wǎng)格變形方法能消除過(guò)量現(xiàn)象，還保持較陡的濃度鋒面。

對(duì)式

設(shè)坐標(biāo)變換

得，

取空間步長(zhǎng)，則

66四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法動(dòng)坐標(biāo)系方法與網(wǎng)格變形方法155156四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法動(dòng)坐標(biāo)系方法與網(wǎng)格變形方法

在固定坐標(biāo)系中取格距為，有

動(dòng)坐標(biāo)系中格點(diǎn)運(yùn)動(dòng)與固定坐標(biāo)系格點(diǎn)運(yùn)動(dòng)重合。

設(shè)固定坐標(biāo)系中格點(diǎn)xj在tn+1時(shí)刻與動(dòng)坐標(biāo)系中Xi重合，有

動(dòng)坐標(biāo)系下濃度結(jié)點(diǎn)的表達(dá)式67四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法動(dòng)坐標(biāo)系方法與網(wǎng)格變形方法156157四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法動(dòng)坐標(biāo)系方法與網(wǎng)格變形方法

式變成

對(duì)半無(wú)限長(zhǎng)砂柱，當(dāng)DL=0時(shí)68四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法動(dòng)坐標(biāo)系方法與網(wǎng)格變形方法157158四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法網(wǎng)格變形方法

設(shè)tn時(shí)刻有限元結(jié)點(diǎn)在空間中的分布是對(duì)應(yīng)濃度分布為

69四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法網(wǎng)格變形方法158159四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法網(wǎng)格變形方法

主要步驟：結(jié)點(diǎn)位置隨時(shí)間變化，基函數(shù)及有限元方程的系數(shù)都依賴于時(shí)間，每推進(jìn)一個(gè)Δt都需要重新計(jì)算。二維問(wèn)題還須進(jìn)行畸形判斷

70四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法網(wǎng)格變形方法159160四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法結(jié)合動(dòng)坐標(biāo)的網(wǎng)格變形方法

一維對(duì)流-彌散在某個(gè)動(dòng)點(diǎn)在t時(shí)刻位置x可寫成

其中，x0是動(dòng)點(diǎn)的初始位置，有取時(shí)間的一階偏導(dǎo)數(shù)對(duì)流-彌散方程寫成

若控制動(dòng)點(diǎn)速度dx/dt并讓它接近流速u，則可轉(zhuǎn)換成彌散為主的的方程式71四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法結(jié)合動(dòng)坐標(biāo)的網(wǎng)格變形方法若160161四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法若使用伽遼金法，基函數(shù)因依賴于結(jié)點(diǎn)位置而成為時(shí)間的函數(shù)。若用表示隨時(shí)間變化的基函數(shù)，則對(duì)任意t，有

故72四、對(duì)流占主導(dǎo)地位時(shí)的數(shù)值方法若使用伽遼金法，基函數(shù)因依161162四、對(duì)流占主