溶質(zhì)運(yùn)移理論(五)水動力彌散方程的數(shù)值解法課件

1、地下水流數(shù)值模擬方法第六章 水動力彌散方程的數(shù)值解法中國地質(zhì)大學(xué)環(huán)境學(xué)院2019春2一、有限差分法基本思想： 將空間劃分成若干網(wǎng)格，把時間分成若干小時段，每一個網(wǎng)格中心點(diǎn)（格點(diǎn)）處的未知變量（H或C）視為該網(wǎng)格平均值。利用差商代替微商?；静襟E： （1）剖分（2）建立差分方程組（3）求解3一、有限差分法-導(dǎo)數(shù)的有限差分近似圖中，去x軸上任意一點(diǎn)i，其坐標(biāo)為 在改點(diǎn)左右相聚為 處分別?。╥-1）和（i+1），其坐標(biāo)分別為 和以i為中心，泰勒展開C(x)4一、有限差分法-導(dǎo)數(shù)的有限差分近似整理并略去余項(xiàng)（6-1）-（6-2），再除以 略去余項(xiàng)（6-1）-（6-2），再除以 略去余項(xiàng)：分別一階導(dǎo)數(shù)的

2、向前差分、向后差分及中心差分二階中心差分5一、有限差分法-一維水動力彌散的差分解法一維水動力彌散方程 （6-7）（1）顯格式式中僅有一個未知數(shù)，解得 式（6-7）中的對流項(xiàng)取中心差分 6可以證明，穩(wěn)定性準(zhǔn)則要求 即 即 。（1）（2）。由條件（2），格距要求很?。挥桑?）可知，鑒于較小，導(dǎo)致不能太大，將（2）代入（1）式中，得到 7若對流項(xiàng)改為向后差分 解得 穩(wěn)定性要求不難看出，穩(wěn)定性限制比對流項(xiàng)取中心差分有所改善。8（2）隱式格式 整理后得到 隱格式是無條件穩(wěn)定的。9（2）Grank-Nicolson格式 整理后得到 取隱式和顯示的平均，即Grank-Nicolson格式10一、有限差分法-

3、二維水動力彌散的差分解法二維水動力彌散方程 （4-56）（1）顯格式式中僅有一個未知數(shù)式（4-56）中的對流項(xiàng)取中心差分 11一、有限差分法-二維水動力彌散的差分解法化簡后，有涉及以（i，j）為中心的5個網(wǎng)格點(diǎn)在tn時刻的已知濃度12一、有限差分法-二維水動力彌散的差分解法（2）隱格式式（4-56）中右端的對流項(xiàng)取中心差分，右端個C的時階均取n+1水平 13一、有限差分法-二維水動力彌散的差分解法（2）隱格式整理后收斂且無條件穩(wěn)定涉及以（i，j）為中心的5個網(wǎng)格點(diǎn)在tn+1時刻的未知濃度14一、有限差分法-二維水動力彌散的差分解法（3）Grank-Nicolson格式將隱式格式的兩式相加除以2

4、15一、有限差分法-二維水動力彌散的差分解法（3）Grank-Nicolson格式整理得16一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱式法（ADI法）分兩次對三對角矩陣求逆，將一個t分成兩個t/2計算 第一個半時間步，對x方向的偏導(dǎo)數(shù)采用隱式差分，對y方向的偏導(dǎo)數(shù)采用顯示差分。17一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱式法（ADI法）整理得 第二個半時間步，對y方向的偏導(dǎo)數(shù)采用隱式差分，對x方向的偏導(dǎo)數(shù)采用顯示差分。18一、有限差分法-二維水動力彌散的差分解法（4）交替方向隱式法（ADI法）整理得收斂且無條件穩(wěn)定19一、有限差分法-求解差分方程的計算機(jī)程序舉例算

5、例 見教材P81-8520二、有限單元法-伽遼金有限單元法伽遼金法 對均勻多孔介質(zhì)，一維穩(wěn)定流場中二維水動力彌散方程 取x軸方向與地下水流向相同，記伽遼金法是尋找一個級數(shù)形式的試探函數(shù)作微分方程的近似解，并使其滿足邊界條件 （6-26）21二、有限單元法-伽遼金有限單元法上式一般不會滿足方程，因?yàn)?僅是近似解，得到一個剩余誤差函數(shù) R(x，y) ，在平均意義下迫使誤差為0 我們加以權(quán)，令剩余的加權(quán)積分為0，W是一組權(quán)函數(shù)。加權(quán)剩余法（6-29）22二、有限單元法-伽遼金有限單元法在整個研究區(qū)D上，基函數(shù)NL(x,y)分段定義(1)函數(shù)區(qū)域連續(xù)（2）一階導(dǎo)數(shù)不連續(xù)（3）二階導(dǎo)不容易確定故采用分步

6、積分的方法使二階導(dǎo)數(shù)降階23二、有限單元法-伽遼金有限單元法對一維積分整理后分步積分推廣到二維的情況（6-32）代入（6-29）24二、有限單元法-伽遼金有限單元法有限單元剖分與基函數(shù)（1）三角形單元 見地下水流動問題數(shù)值模擬（2）矩形單元將區(qū)域記為D，邊界記為B。要求各單元均質(zhì)、等厚，即T、為常數(shù)。結(jié)點(diǎn)（內(nèi)結(jié)點(diǎn)、邊界結(jié)點(diǎn)） 25二、有限單元法-伽遼金有限單元法構(gòu)造單元函數(shù)NL基函數(shù)取為 “雙線性插值” 基函數(shù)NL（x，y）形狀如同一頂高等于1、有4條直線斜邊和4條下凹型曲邊的尖頂斗笠26二、有限單元法-伽遼金有限單元法子區(qū)DL以結(jié)點(diǎn)L為其共同結(jié)點(diǎn)的所有矩形單元（0）27二、有限單元法-伽遼金

7、有限單元法典型矩形單元 該單元中心位于坐標(biāo)原點(diǎn)處，且4條邊分別平行x，y軸 ,結(jié)點(diǎn)從左下角開始按逆時針方向編號i,j,m,n。沿x方向長2a，y方向?qū)?b。28二、有限單元法-伽遼金有限單元法按上述要求所構(gòu)造的NeL形式為對典型單元，令 則29二、有限單元法-伽遼金有限單元法其中即（6-34）30二、有限單元法-伽遼金有限單元法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)濃度用 表示31二、有限單元法-伽遼金有限單元法結(jié)點(diǎn)濃度CL(t)和單元基函數(shù)NeL(x,y)來確定單元內(nèi)的近似解 ，即對

8、于區(qū)域D，結(jié)點(diǎn)L的基函數(shù)在子區(qū)DL內(nèi)各單元分塊確定。稱為矩陣單元e上的基函數(shù)（6-36）32二、有限單元法-伽遼金有限單元法伽遼金有限單元法 左端積分范圍為D，由于在 上NL=0,改寫在子區(qū)DL的積分，對矩形單元逐個積分、求和。 設(shè)子區(qū)DL上有neL個單元 ，有（6-39）33二、有限單元法-伽遼金有限單元法由（6-36）得由（6-34）得再由（6-36）得（6-40）（6-42）34二、有限單元法-伽遼金有限單元法將（6-40） （6-42） 代入（6-39）得35二、有限單元法-伽遼金有限單元法（6-44）36二、有限單元法-伽遼金有限單元法用矩陣形式表示一階微分方程組將一階導(dǎo)數(shù)用差分形式

9、表示按矩陣符號有若C在新時間水平t+t上取值，則37二、有限單元法-伽遼金有限單元法（1）給定初始條件和邊界條件，可求（6-46）各個系數(shù)矩陣，于是可解得第一個 末時刻各結(jié)點(diǎn)的濃度 ；（2）以此濃度為已知濃度 ，計算新的系數(shù)矩陣 ；（3）再次求解，得到第二個時間步長末時刻結(jié)點(diǎn)濃度（6-46）濃度 可在 之間任意位置近似表示38二、有限單元法-伽遼金有限單元法按單元順序計算系數(shù)矩陣每個矩陣單元涉及4個結(jié)點(diǎn)i，j，m，n，涉及4個方程單元e的作用，可用4行與4列的矩陣來表示，稱為單元矩陣39二、有限單元法-伽遼金有限單元法算例：單元彌散矩陣與整體彌散矩陣的關(guān)系設(shè)研究區(qū)劃分為8個單元，16個結(jié)點(diǎn)，四

10、周為隔水邊界編號如下： 40二、有限單元法-伽遼金有限單元法算例：單元彌散矩陣與整體彌散矩陣的關(guān)系計算每個單元彌散矩陣，按雙下標(biāo)已知，放置單元彌散矩陣并疊加，得到總體彌散矩陣。41二、有限單元法-伽遼金有限單元法從（6-44）知將坐標(biāo)原點(diǎn)移到矩形單元中心有42二、有限單元法-伽遼金有限單元法通過雙下標(biāo)一致的元素求和，得到總體彌散矩陣諸元素43二、有限單元法-伽遼金有限單元法積分求解可用高斯求積方法，詳見P97-10544二、有限單元法-伽遼金有限單元法45二、有限單元法-伽遼金有限單元法算例：二維水動力彌散的伽遼金法計算機(jī)程序 數(shù)學(xué)模型寫成剖分成9個單元共20個結(jié)點(diǎn)，空間步長取x=y=10m，

11、時間步長t=5d，滲透流速u=0.01m/d 具體代碼見P106-11146三、過程與數(shù)值彌散過量現(xiàn)象 如一維流動水動力彌散中單位濃度梯度函數(shù)注入的情況。 對比解析解與數(shù)值解C-x曲線，在濃度封面附近數(shù)值計算的濃度超過最大濃度值1和小于最小濃度0，失去物理意義。數(shù)值彌散 若純對流方程采用有限差分逼近，其結(jié)果呈現(xiàn)似有彌散的存在。47三、過程與數(shù)值彌散48三、過程與數(shù)值彌散過量現(xiàn)象-解析 由于時間步長與空間尺度匹配不佳，因而含水層在數(shù)值上不能“吸收”住污染物所引起。解決辦法： 將時間增量選擇為幾何級數(shù)的子項(xiàng)； 做試驗(yàn)校正時間步長和差分網(wǎng)格；49三、過程與數(shù)值彌散數(shù)值彌散-解析 由截斷誤差引起。 對

12、流項(xiàng)一階空間導(dǎo)數(shù)采用向后差分逼近，濃度Ci-1泰勒展開 變形得 向后差分的截斷誤差主部分為相當(dāng)于彌散系數(shù)的彌散項(xiàng)50三、過程與數(shù)值彌散數(shù)值彌散-解析 由故純對流方程有限差分近似表示后，結(jié)果如同對流-彌散方程的結(jié)果，記用數(shù)值彌散系數(shù)與物理彌散系數(shù)的比值來評價數(shù)值彌散對物理彌散的干涉程度此條件下的數(shù)值彌散系數(shù)51三、過程與數(shù)值彌散數(shù)值彌散-解析 定義為Peclet數(shù)， 時，彌散占優(yōu)； 時，對流占優(yōu)。速度u與物理彌散系數(shù)DL比越大， 越大。對于小 數(shù)，數(shù)值解與精確解非常接近，對于大 數(shù)，數(shù)值解過渡帶變寬，濃度鋒面變平緩，在鋒面的前后還出現(xiàn)解的振動過量現(xiàn)象。控制x的大小從而減少數(shù)值彌散52三、過程與數(shù)

13、值彌散對純對流方程 53三、過程與數(shù)值彌散對純對流方程 結(jié)果也造成數(shù)值彌散此時，總的數(shù)值彌散是兩者之和54四、對流占主導(dǎo)地位時的數(shù)值方法上游加權(quán)法 對流項(xiàng)的有限差分化過程中乘上一個適當(dāng)?shù)臋?quán)因子，則波動和過量得到改善或消除一維水動力彌散中對流項(xiàng)取中心隱式差分若相鄰格點(diǎn)i-1或（和）i+1的濃度越大，則格點(diǎn)i的濃度也越大，反映到差分方程上，系數(shù)Gi-1和Gi+1不得與Gi同號，上式要求 不能滿足時，對流占優(yōu)勢55四、對流占主導(dǎo)地位時的數(shù)值方法上游加權(quán)法 假定DL和u為常數(shù)，取等格距差分網(wǎng)格。一維對流-彌散方程在 積分，有對流-彌散通量用有限差分近似表示 56四、對流占主導(dǎo)地位時的數(shù)值方法上游加權(quán)法

14、 為權(quán)因子；當(dāng)=1/2時，對流項(xiàng)相當(dāng)于中心差分；當(dāng)=1時，對流項(xiàng)相當(dāng)于取向后差分， 當(dāng)=0時，對流項(xiàng)相當(dāng)于向前差分。速度不同， 取值不同57四、對流占主導(dǎo)地位時的數(shù)值方法上游加權(quán)法 上游濃度權(quán)因子大于1/2。物理意義：加強(qiáng)上游格點(diǎn)濃度值在對流項(xiàng)中的作用可通過濃度的線性（或二次）插值具體確定在t期間C在格邊上隨時間的變化，從而確定值。58四、對流占主導(dǎo)地位時的數(shù)值方法 仍用Gi-1、Gi和Gi+1分別表示上式左端第一、二、三項(xiàng)括號內(nèi)系數(shù)。則采用上游加權(quán)法有限差分格式，可消除過量和波動，但截斷誤差增大，數(shù)值彌散反而增大。中心差分格式，截斷誤差為二階上游加權(quán)有限差分截斷誤差為一階59四、對流占主導(dǎo)地

15、位時的數(shù)值方法 上游權(quán)因子的物理實(shí)質(zhì)：對流項(xiàng)濃度C在某斷面，例如 斷面的中心差分意味著對流項(xiàng)的濃度在t時段內(nèi)均以i與i-1兩格點(diǎn)的平均濃度，即 通過 斷面；實(shí)際為： t時段的初始時刻是該平均濃度值，隨后在對流的作用下，該斷面的濃度顯然會向上游格點(diǎn)的濃度Ci-1 (u0時)的變化。60四、對流占主導(dǎo)地位時的數(shù)值方法特征值方法 令源匯項(xiàng)為0，展開對流-彌散方程 61四、對流占主導(dǎo)地位時的數(shù)值方法特征值方法 假定流體密度值不變（不可壓縮），含水層不可壓縮 ，則 ，（6-75）簡化為基本思路：物質(zhì)導(dǎo)數(shù)表示一個沿軌跡（特征曲線） 運(yùn)動的研究對象濃度變化率。由彌散作用引起。若存在運(yùn)會想，還將受到沿軌跡源與

16、匯、化學(xué)、生物化學(xué)反應(yīng)、吸附與解析等。對流項(xiàng)用沿特征線的示蹤質(zhì)點(diǎn)的運(yùn)動描述 （6-76）62四、對流占主導(dǎo)地位時的數(shù)值方法特征值方法 特征點(diǎn)均勻分布研究區(qū)，賦一個初始濃度。 含水層中未受污染部分特征點(diǎn)的初始濃度為零，特征點(diǎn)濃度眼軌跡的瞬時濃度變化，借助矩形網(wǎng)格差分格式計算。63四、對流占主導(dǎo)地位時的數(shù)值方法64四、對流占主導(dǎo)地位時的數(shù)值方法模型精確度 取決于網(wǎng)格距大小?？梢栽谕痪W(wǎng)格上計算速度場。如果采用較少的網(wǎng)格距，可以改善流場的描述。 坐標(biāo)為 的特征點(diǎn)所在的單元，可用下標(biāo)i,j確定 精確度受到每一網(wǎng)格內(nèi)所選擇的最初的特征點(diǎn)數(shù)目影響。65四、對流占主導(dǎo)地位時的數(shù)值方法缺點(diǎn)： （1）實(shí)際程序冗

17、長，需記錄特征點(diǎn)的蹤跡； （2）抽水點(diǎn)或出流邊界的特征點(diǎn)需要消除，入滲邊界須有特征點(diǎn)生成，不透水邊界特征點(diǎn)須反射回來； （3）要通過消去特征點(diǎn)來避免滯留點(diǎn)附近特征點(diǎn)堆積； （4）特征點(diǎn)耗盡的單元，必須產(chǎn)生新的特征點(diǎn)。66四、對流占主導(dǎo)地位時的數(shù)值方法動坐標(biāo)系方法與網(wǎng)格變形方法 能消除過量現(xiàn)象，還保持較陡的濃度鋒面。 對式 設(shè)坐標(biāo)變換 得， 取空間步長 ，則 67四、對流占主導(dǎo)地位時的數(shù)值方法動坐標(biāo)系方法與網(wǎng)格變形方法 在固定坐標(biāo)系中取格距為 ，有 動坐標(biāo)系中格點(diǎn)運(yùn)動與固定坐標(biāo)系格點(diǎn)運(yùn)動重合。 設(shè)固定坐標(biāo)系中格點(diǎn)xj在tn+1時刻與動坐標(biāo)系中Xi重合，有 動坐標(biāo)系下濃度結(jié)點(diǎn)的表達(dá)式68四、對流占

18、主導(dǎo)地位時的數(shù)值方法動坐標(biāo)系方法與網(wǎng)格變形方法 式變成 對半無限長砂柱，當(dāng)DL=0時69四、對流占主導(dǎo)地位時的數(shù)值方法網(wǎng)格變形方法 設(shè)tn時刻有限元結(jié)點(diǎn)在空間中的分布是對應(yīng)濃度分布為 70四、對流占主導(dǎo)地位時的數(shù)值方法網(wǎng)格變形方法 主要步驟：結(jié)點(diǎn)位置隨時間變化，基函數(shù)及有限元方程的系數(shù)都依賴于時間，每推進(jìn)一個 t都需要重新計算。二維問題還須進(jìn)行畸形判斷 71四、對流占主導(dǎo)地位時的數(shù)值方法結(jié)合動坐標(biāo)的網(wǎng)格變形方法 一維對流-彌散在某個動點(diǎn)在t時刻位置x可寫成 其中，x0是動點(diǎn)的初始位置，有取時間的一階偏導(dǎo)數(shù)對流-彌散方程寫成 若控制動點(diǎn)速度dx/dt 并讓它接近流速u，則可轉(zhuǎn)換成彌散為主的的方程

19、式72四、對流占主導(dǎo)地位時的數(shù)值方法若使用伽遼金法，基函數(shù)因依賴于結(jié)點(diǎn)位置而成為時間的函數(shù)。若用 表示隨時間變化的基函數(shù)，則對任意t，有 故73四、對流占主導(dǎo)地位時的數(shù)值方法隨機(jī)步行法 采用示蹤劑描述污染物運(yùn)移。只有被污染的特征點(diǎn)在流場中運(yùn)動。每個特征點(diǎn)都給定一個不變的污染物質(zhì)量，此質(zhì)量的和等于排進(jìn)含水層的總污染物質(zhì)量。 取許多單個特征點(diǎn)軌跡（隨機(jī)步行）的平均值，可得到一個有彌散的特征點(diǎn)的分布。要得到一個濃度分布，就要先疊加一個網(wǎng)格并算每個網(wǎng)格單元所含有的特征點(diǎn)數(shù)。 74四、對流占主導(dǎo)地位時的數(shù)值方法隨機(jī)步行法-以一維運(yùn)移為例在x=0處，一個瞬時注入的、質(zhì)量為 M的理想失蹤及，在t時刻濃度分布

20、對一個固定時間t，看做正態(tài)分布 75四、對流占主導(dǎo)地位時的數(shù)值方法隨機(jī)步行法-以一維運(yùn)移為例在時間t=0和位置x=0處，放置具有污染物質(zhì)量為M/N的N個特征點(diǎn)，每一個特征點(diǎn)移動距離x，而在時刻t到達(dá)它們各自位置，即Z是一個平均值為0，標(biāo)準(zhǔn)差為1的正態(tài)分布隨機(jī)變量 所得的頻率分布，通過一個標(biāo)準(zhǔn)化因子，有對有限數(shù)量的顆粒，可近似去頂C(x,t)，將空間變量劃分為長度x的若干區(qū)間，區(qū)間中點(diǎn)濃度76四、對流占主導(dǎo)地位時的數(shù)值方法隨機(jī)步行法-以一維運(yùn)移為例若孔隙平均流速是時間和位置的函數(shù)，那么時間為0開始的顆粒P在時間間隔t內(nèi)的隨機(jī)軌跡xp(t)為 77四、對流占主導(dǎo)地位時的數(shù)值方法隨機(jī)步行法-二維流場

21、 首先假定流動方向平行于x軸，得到Z和Z是正態(tài)分布的隨機(jī)變量的兩個值。對任意方向的流動，需考慮彌散張量特性。78四、對流占主導(dǎo)地位時的數(shù)值方法隨機(jī)步行法-二維流場 首先假定流動方向平行于x軸，得到Z和Z是正態(tài)分布的隨機(jī)變量的兩個值。對任意方向的流動，需考慮彌散張量特性。對流運(yùn)動可確定在水流流動方向及其正交方向的彌散運(yùn)動79四、對流占主導(dǎo)地位時的數(shù)值方法隨機(jī)步行法 80四、對流占主導(dǎo)地位時的數(shù)值方法隨機(jī)步行法 數(shù)值結(jié)果的質(zhì)量主要取決于特征點(diǎn)數(shù)目N，要考慮網(wǎng)格的選擇。大網(wǎng)格間距會產(chǎn)生極平均的結(jié)果，小網(wǎng)格間距會產(chǎn)生十分粗糙的分布。 在穩(wěn)定流動和一個源或數(shù)個源都有固定強(qiáng)度的情況下，可以節(jié)省很多特征點(diǎn)。8182四、對流占主導(dǎo)地位時的數(shù)值方法隨機(jī)步行法