陳榮三， 王凡芮， 鄒 敏

(中國(guó)地質(zhì)大學(xué)(武漢) 數(shù)學(xué)與物理學(xué)院， 湖北 武漢 430074)

對(duì)流-彌散方程是地下水溶質(zhì)運(yùn)移的控制方程.簡(jiǎn)單的理想的對(duì)流-彌散方程可以求出解析解，對(duì)于復(fù)雜的方程還要借助于數(shù)值方法.近年來，國(guó)內(nèi)外學(xué)者針對(duì)對(duì)流-彌散方程提出了大量的數(shù)值方法，大致可以分成三類：Euler法、Lagrange法和混合Euler-Lagrange法.這些方法都能較好地解決彌散占優(yōu)的問題.但是一旦對(duì)流占優(yōu)，數(shù)值解會(huì)出現(xiàn)過量問題、過多的數(shù)值彌散和非物理振蕩問題.

有限差分法、有限體積法和有限元法是三種比較常見的Euler方法.Roache[1]提出了上游加權(quán)法，該方法可以消除非物理振蕩，但是增大了數(shù)值彌散.隨后許多學(xué)者[2-3]提出了不同的上游加權(quán)法，上游加權(quán)對(duì)減小解的振蕩作用有限, 但卻會(huì)減低解的精度.為了提高解的精度,Brooks等[4]提出了沿流線上游加權(quán)的彼得羅夫-迦遼金方法.彼得羅夫-迦遼金方法在消除了數(shù)值振動(dòng)現(xiàn)象的同時(shí)增大了數(shù)值彌散.梅一等[5]發(fā)展了一種帶彌散因子的上游加權(quán)法，通過調(diào)整彌散因子可以控制非物理振蕩但不降低解的精度.Leisman等[6]提出了引入人工擴(kuò)散量和加權(quán)法，可以使得系數(shù)矩陣對(duì)稱正定.劉揚(yáng)[7]發(fā)展了帶加罰函數(shù)的改進(jìn)的Crank-Nicholson方法，該方法可以保證解的穩(wěn)定性同時(shí)不會(huì)出現(xiàn)非物理振蕩.成建梅等[8]總結(jié)了飽和水流中溶質(zhì)運(yùn)移方程求解的各種數(shù)值方法.焦甜等[9]分析了Crank-Nicholson方法，上游加權(quán)法和人工增量法的截?cái)嗾`差、穩(wěn)定性和收斂性.Cheng等[10]將局部間斷有限元方法推廣到了對(duì)流-彌散方程的計(jì)算.

近幾年來，茅德康等[11-16]針對(duì)雙曲守恒性方程提出了一種Godunov型的熵格式.與傳統(tǒng)Godunov不同的是，熵格式除了數(shù)值計(jì)算數(shù)值解外，還數(shù)值計(jì)算數(shù)值熵，通過數(shù)值模擬多個(gè)物理量可以從多個(gè)側(cè)面來刻畫物理過程.數(shù)值實(shí)驗(yàn)證明該方法可以有效地控制數(shù)值黏性并且適合長(zhǎng)時(shí)間的數(shù)值計(jì)算，在保證數(shù)值穩(wěn)定的同時(shí)最大程度上減少了數(shù)值彌散.崔艷芬等[10]對(duì)熵格式進(jìn)行了數(shù)值分析，揭示了這種格式具有在計(jì)算中各步誤差相互抵消的性質(zhì).

本文將熵格式推廣到地下水溶質(zhì)運(yùn)移方程.首先采用分裂方法將地下水溶質(zhì)運(yùn)移方程分成對(duì)流方 程和彌散方程，對(duì)流方程是一個(gè)雙曲型方程，利用熵格式求解，彌散方程的空間離散用二階中心格式離散,時(shí)間離散用簡(jiǎn)單的向前差分離散.通過數(shù)值實(shí)驗(yàn)，對(duì)不同對(duì)流強(qiáng)度的地下水溶質(zhì)運(yùn)移方程進(jìn)行了數(shù)值計(jì)算，計(jì)算結(jié)果表明,熵格式?jīng)]有出現(xiàn)過量問題，沒有出現(xiàn)非物理振蕩，數(shù)值彌散小，特別適合強(qiáng)對(duì)流問題的數(shù)值計(jì)算.

1 熵格式的描述

考慮如下的地下水溶質(zhì)運(yùn)移方程[9]

(1)

式中:ρ為質(zhì)量濃度；DL為水動(dòng)力彌散系數(shù)；u為速度；x為距離；t為時(shí)間.在式(1)兩邊同時(shí)乘以2ρ,整理得

(2)

這里U(ρ)=ρ2，把它稱之為熵.熵格式是Godunov型有限體積格式，與傳統(tǒng)的Godunov格式不同的是，熵格式有兩個(gè)數(shù)值實(shí)體，即數(shù)值解和數(shù)值熵.也就是說，不僅要計(jì)算水溶質(zhì)運(yùn)移方程(1)還需要計(jì)算相應(yīng)的熵方程(2).為了描述的簡(jiǎn)單，采用一致等分網(wǎng)格，網(wǎng)格的寬度為Δx，時(shí)間步長(zhǎng)為Δt，第i網(wǎng)格單元為Ii=(xi-1/2,xi+1/2)，其中xi±1/2=xi±0.5Δx，xi表示網(wǎng)格中心.引入數(shù)值解ρi，n為網(wǎng)格單元Ii在tn時(shí)刻精確解的網(wǎng)格平均的近似，即

(3)

另外引入了一個(gè)數(shù)值熵Ui,n，它是對(duì)精確解的熵的網(wǎng)格平均的逼近，即

(4)

采用分裂的方法來求解式(1)和式(2)，將式(1)和式(2)分成對(duì)流部分和彌散部分.式(1)和式(2)的對(duì)流部分和彌散部分分別為

(5)

和

(6)

利用熵格式計(jì)算對(duì)流方程(5)，熵格式包含數(shù)值解ρi,n和數(shù)值熵Ui,n兩個(gè)數(shù)值實(shí)體. 重構(gòu)函數(shù)和ρi,n和Ui,n都有關(guān).熵格式按照重構(gòu)、發(fā)展和網(wǎng)格平均三步來進(jìn)行.具體步驟如下：

(1) 重構(gòu).tn時(shí)刻第i個(gè)網(wǎng)格Ii內(nèi)采用熵臺(tái)階重構(gòu)，它是一個(gè)階梯函數(shù)，如式(7)所示.

(7)

其中di,n為重構(gòu)斜率. 顯然，熵臺(tái)階重構(gòu)式(7)自動(dòng)滿足

(8)

首先計(jì)算熵臺(tái)階di,n,e，要求

(9)

式(9)意味著重構(gòu)解的熵的網(wǎng)格平均等于該網(wǎng)格的數(shù)值熵. 將式(7)代入式(9)可以計(jì)算出熵臺(tái)階di,n,e

(10)

di,n,e符號(hào)取為和(ρj+1,n-ρj-1,n)相同. 另外，定義一個(gè)總變差下降(TVD)的臺(tái)階di,n,TVD

di,n,TVD=min mod(ρi,n-ρi-1,n,ρi+1,n-ρi,n)

(11)

這里min mod定義為

min mod(a,b)=

(12)

最終的臺(tái)階取為di,n,e和di,n,TVD的最小值，符號(hào)取為和(ρi+1,n-ρi-1,n)相同，即

di,n=sgn(ρi+1-ρi-1)min(|di,n,e|,|di,n,TVD|)

(13)

(2) 發(fā)展.以重構(gòu)函數(shù)R(x;ρn,Un)作為tn層的初值，求解以下初值問題:

(14)

定解問題(14)的精確解為v(x,t)=R(x-ut;ρn,Un).

(3) 網(wǎng)格平均. 在t=tn+1時(shí)的數(shù)值解和數(shù)值熵分別計(jì)算為

(15)

(16)

利用中心格式計(jì)算彌散方程(6)中的二階導(dǎo)數(shù)，利用Euler向前差分進(jìn)行時(shí)間離散.最終的格式為

(17)

(ρi+1,n-2ρi,n+ρi-1,n)

(18)

2 數(shù)值試驗(yàn)和結(jié)果分析

算例1 均勻流場(chǎng)[9].考慮如下地下水溶質(zhì)運(yùn)移問題，假設(shè)半無界均質(zhì)砂礫，地下水的速度u是常數(shù)，初始質(zhì)量濃度為0,左邊界為質(zhì)量濃度邊界，右邊界為0邊界，即右無窮遠(yuǎn)處的質(zhì)量濃度為0.該問題的解析解為

(19)

式中:ρ0=1 mg·L-1;u=6 m·d-1;DL=uα，α是縱向彌散率.

Peclet數(shù)定義為Pe=Δx/α，Pe是用來描述對(duì)流和彌散的相對(duì)大小，根據(jù)Pe的大小可以將對(duì)流-彌散方程分為兩大類：當(dāng)Pe>2時(shí)，稱之為對(duì)流占優(yōu)方程；當(dāng)Pe≤2 時(shí)，稱之為彌散占優(yōu)方程.計(jì)算區(qū)域?yàn)閇0,60] m，網(wǎng)格數(shù)取為30.圖1～ 4是對(duì)流強(qiáng)度Pe分別為1、4、16、32，時(shí)間分別為1 d和4 d的數(shù)值解和解析解的質(zhì)量濃度分布對(duì)比圖，圖中圓圈表示數(shù)值解，實(shí)線表示解析解.從圖1～4中不難看出，熵格式?jīng)]有出現(xiàn)過量問題，也沒有出現(xiàn)非物理振蕩，所產(chǎn)生的數(shù)值彌散比較小.隨著對(duì)流強(qiáng)度的增加，數(shù)值解和解析解吻合得更好，這說明熵格式在處理強(qiáng)對(duì)流問題方面具有很大優(yōu)勢(shì).

圖1 Pe=1時(shí)質(zhì)量濃度分布

算例2 非均勻流場(chǎng)[7].考慮如下對(duì)流擴(kuò)散問題

(20)

該問題的精確解為x(1-x)et.圖5是t=0.6的數(shù)值解和解析解的對(duì)比圖，圖中圓圈表示數(shù)值解，實(shí)線表示解析解.從圖5不難看出，數(shù)值解和精確解幾乎重合，沒有出現(xiàn)非物理振蕩也沒出現(xiàn)過量現(xiàn)象.

圖5 t=0.6時(shí)質(zhì)量濃度分布

3 結(jié)論

地下水運(yùn)移方程可以用對(duì)流-彌散方程來刻畫，當(dāng)對(duì)流占優(yōu)的時(shí)候，數(shù)值解容易出現(xiàn)過量和非物理振蕩，一些能克服非物理振蕩的格式往往會(huì)引入過多的數(shù)值彌散，使得解的精度降低.本文將針對(duì)雙曲守恒型方程的熵格式推廣到地下水運(yùn)移方程的計(jì)算，具體做法是采用分裂的方法將地下水運(yùn)移方程分成對(duì)流方程和彌散方程，對(duì)流方程利用熵格式求解，彌散方程用中心格式進(jìn)行求解.數(shù)值試驗(yàn)表明,本文的格式不會(huì)產(chǎn)生過量問題也不會(huì)出現(xiàn)非物理振蕩，特別適合強(qiáng)對(duì)流問題的數(shù)值計(jì)算.