地下水溶質(zhì)運(yùn)移方程有限差分格式的實(shí)證研究

焦　甜,王軍霞,唐仲華，劉　耿

(1.中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院，湖北 武漢 430074；2.中國(guó)地質(zhì)大學(xué)(武漢)數(shù)學(xué)與物理學(xué)院，湖北 武漢 430074；3.北京得力合環(huán)境治理有限公司，北京 100025)

?

地下水溶質(zhì)運(yùn)移方程有限差分格式的實(shí)證研究

焦甜1,王軍霞2,唐仲華1，劉耿3

(1.中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院，湖北 武漢 430074；2.中國(guó)地質(zhì)大學(xué)(武漢)數(shù)學(xué)與物理學(xué)院，湖北 武漢 430074；3.北京得力合環(huán)境治理有限公司，北京 100025)

摘要：在求解對(duì)流占優(yōu)的地下水溶質(zhì)運(yùn)移方程時(shí)，為克服有限差分法和有限單元法產(chǎn)生的過(guò)量和值彌散等虛假數(shù)值現(xiàn)象，出現(xiàn)了許多數(shù)值方法，其中包括Crank-Nicholson格式、上游加權(quán)法和人工擴(kuò)散量法等。以一維地下水溶質(zhì)運(yùn)移方程為例，分析三種差分格式的截?cái)嗾`差、穩(wěn)定性和收斂性，并通過(guò)數(shù)值試驗(yàn)，在不同的對(duì)流強(qiáng)度下，將不同算法的數(shù)值解與解析解進(jìn)行對(duì)比，研究不同算法的優(yōu)缺點(diǎn)和適用條件。結(jié)果表明：在對(duì)流占優(yōu)較弱時(shí)，三種方法都適用，但Crank-Nicholson差分格式誤差最??；當(dāng)對(duì)流強(qiáng)度增大時(shí)，適用時(shí)間步長(zhǎng)較小的人工擴(kuò)散量法；當(dāng)對(duì)流強(qiáng)度較大時(shí)，適用時(shí)間步長(zhǎng)較大的人工擴(kuò)散量法，或者空間步長(zhǎng)較小的上游加權(quán)法。該研究可為解決實(shí)際問(wèn)題時(shí)算法的選擇提供依據(jù)。

關(guān)鍵詞：地下水溶質(zhì)運(yùn)移方程；對(duì)流占優(yōu)強(qiáng)度；有限差分格式；解析解；實(shí)證對(duì)比

隨著經(jīng)濟(jì)的快速發(fā)展，地下水污染問(wèn)題日益嚴(yán)重，由于地下水污染具有滯后性、隱蔽性以及一旦污染便難以清除的特點(diǎn)[1]，因此地下水污染問(wèn)題得到了廣泛關(guān)注。盧曉華[2]通過(guò)建立地下水溶質(zhì)運(yùn)移模型，對(duì)廢水處理車間發(fā)生泄漏情況下地下水重金屬污染情景進(jìn)行了數(shù)值模擬與環(huán)境影響預(yù)測(cè)評(píng)價(jià)，對(duì)環(huán)保部門(mén)從源頭上防治和管理潛在的地下水重金屬污染、保護(hù)地下水環(huán)境具有重要的意義；李合蓮等[3]采用有限元法對(duì)尾礦中238U在地下水中的遷移情況進(jìn)行了模擬；鄧鼎興[4]采用有限元法對(duì)礦山地下水污染的分區(qū)治理效果進(jìn)行了數(shù)值模擬檢驗(yàn)，結(jié)果表明數(shù)值模擬預(yù)測(cè)結(jié)果與實(shí)際監(jiān)測(cè)值較為接近，可為類似礦山地下水污染治理提供借鑒。但地下水溶質(zhì)運(yùn)移方程的求解精度會(huì)直接影響對(duì)水環(huán)境的分析，例如當(dāng)對(duì)流占優(yōu)時(shí)，傳統(tǒng)的數(shù)值方法(如有限單元法和有限差分法)會(huì)產(chǎn)生數(shù)值彌散和數(shù)值波動(dòng)等虛假數(shù)值現(xiàn)象。目前，針對(duì)對(duì)流作用占優(yōu)的地下水溶質(zhì)運(yùn)移方程求解已經(jīng)有許多方法，包括歐拉法、拉格朗日法和混合-歐拉拉格朗日法[5]等，但這些方法仍然存在許多科學(xué)問(wèn)題亟待解決。

有限差分法和有限單元法[6]等是常用的歐拉法。對(duì)于有限差分法，有學(xué)者提出了上游加權(quán)的方法,此后不少學(xué)者也對(duì)此進(jìn)行了研究[7]，結(jié)果表明此方法消除了數(shù)值波動(dòng)，但卻增大了數(shù)值彌散。相比之下有限單元法精度高、較靈活，但當(dāng)對(duì)流占優(yōu)時(shí)，時(shí)間導(dǎo)數(shù)若采用中心差分會(huì)出現(xiàn)過(guò)量，向后差分則會(huì)出現(xiàn)數(shù)值彌散，對(duì)此很多學(xué)者提出了不同的上游加權(quán)方法[8]，但由于上游加權(quán)法結(jié)點(diǎn)所在的邊線不一定是流線方向，就出現(xiàn)了“側(cè)風(fēng)現(xiàn)象”，為了消除這種現(xiàn)象，Brooks等[9]提出了沿流線上游加權(quán)的彼得羅夫-迦遼金方法，但實(shí)踐表明這種方法也是在消除數(shù)值振蕩的同時(shí)導(dǎo)致了更大的數(shù)值彌散。梅一等[10]在上游加權(quán)法的基礎(chǔ)上引入一個(gè)數(shù)值彌散因子，并適當(dāng)調(diào)節(jié)因子的大小來(lái)有效控制數(shù)值振蕩，減少數(shù)值彌散,并且達(dá)到了滿意的精度。動(dòng)坐標(biāo)法通過(guò)控制坐標(biāo)系的運(yùn)動(dòng)速度有效減小了求解對(duì)流項(xiàng)引起的數(shù)值誤差;網(wǎng)格變形法通過(guò)固定坐標(biāo)系，使網(wǎng)格結(jié)點(diǎn)自動(dòng)跟著鋒面的移動(dòng)而聚在鋒面附近，來(lái)減小鋒面附近的Peclet值，消除過(guò)量和數(shù)值彌散問(wèn)題。雖然這兩種方法處理對(duì)流占優(yōu)問(wèn)題時(shí)比較有效，并且保持了較高的精度，但實(shí)際問(wèn)題中計(jì)算工作量比較大，難以推廣到三維問(wèn)題上。由于對(duì)流項(xiàng)占主導(dǎo)地位時(shí)，地下水溶質(zhì)運(yùn)移方程接近于一個(gè)雙曲線方程，所以啟發(fā)人們利用雙曲線的特征線法來(lái)求解，Garder首先用此法求解可溶混溶質(zhì)運(yùn)移遷移，此種方法原理簡(jiǎn)單、程序冗長(zhǎng)，但能有效減少數(shù)值彌散，也能處理各種復(fù)雜條件。Douglas等[11]提出了特征有限元法和有限差分法，其數(shù)值計(jì)算方面的有效性曾引起了許多學(xué)者的關(guān)注，并得到了廣泛的應(yīng)用研究[12]。隨機(jī)步行法將溶質(zhì)運(yùn)移視為大量示蹤劑特征點(diǎn)運(yùn)動(dòng)的平均結(jié)果，當(dāng)特征點(diǎn)足夠多時(shí)，基本上能刻畫(huà)溶質(zhì)運(yùn)移過(guò)程，優(yōu)點(diǎn)是由于特征點(diǎn)在研究區(qū)域內(nèi)連續(xù)運(yùn)動(dòng)，所以計(jì)算濃度時(shí)只需要統(tǒng)計(jì)區(qū)域上特征點(diǎn)的數(shù)目，不用在每個(gè)時(shí)間步長(zhǎng)內(nèi)計(jì)算整個(gè)研究區(qū)的濃度。Kinzellbach將隨機(jī)步行法與有限差分法、有限單元法和特征值法進(jìn)行了對(duì)比，開(kāi)發(fā)出與MODFLOW配合使用的MODWALK軟件，但是其計(jì)算精度取決于示蹤點(diǎn)數(shù)目，如果數(shù)目不夠，統(tǒng)計(jì)出的濃度解會(huì)比較粗糙。Leismann等[13]提出了引入人工擴(kuò)散量和加權(quán)法來(lái)求解，定義并引入人工擴(kuò)散項(xiàng)進(jìn)行校正。劉揚(yáng)[14]提出一種加罰函數(shù)的Crank-Nicholson差分格式，形成的代數(shù)矩陣對(duì)角元嚴(yán)格占優(yōu)，在保持?jǐn)?shù)值解穩(wěn)定性的基礎(chǔ)上有效克服了數(shù)值振蕩。

對(duì)于各種方法的優(yōu)劣性，成建梅等[15]對(duì)溶質(zhì)運(yùn)移問(wèn)題數(shù)值解從歐拉法、拉格朗日法、混合-歐拉拉格朗日法三大方面進(jìn)行了現(xiàn)狀綜述；劉苑等[16]對(duì)比了MT3D模塊中MOC、MMOC、HMOC和UFDM 4種數(shù)值方法；寧一鳴[17]對(duì)比了二維問(wèn)題的彌散占優(yōu)和對(duì)流占優(yōu)的幾種數(shù)值解差分格式，但沒(méi)有討論不同對(duì)流強(qiáng)度的情況，并探討了一種基于三角形網(wǎng)格的上游加權(quán)有限差分格式。除此之外，許多學(xué)者也從各個(gè)方面研究了對(duì)流占優(yōu)數(shù)值解差分格式的特點(diǎn)，但是對(duì)不同對(duì)流占優(yōu)強(qiáng)度下不同數(shù)值解差分格式的適用性沒(méi)有做詳細(xì)的討論。因此，本文在前人研究工作的基礎(chǔ)上，對(duì)比了三種差分格式的截?cái)嗾`差、穩(wěn)定性和收斂性，并通過(guò)數(shù)值試驗(yàn)在不同的對(duì)流強(qiáng)度下，將不同算法的數(shù)值解與解析解進(jìn)行對(duì)比，研究了不同算法的優(yōu)缺點(diǎn)和適用條件，討論了時(shí)間步長(zhǎng)對(duì)解的影響，以便于研究者在實(shí)際問(wèn)題中選擇適合的算法以及網(wǎng)格剖分大小，從而提高模擬的精度。

1三種差分格式

對(duì)于一維對(duì)流-彌散問(wèn)題[18]，地下水溶質(zhì)運(yùn)移方程為

(1)

式中：C為濃度;DL為水動(dòng)力彌散系數(shù);u為地下水流速;x為空間距離;t為時(shí)間。

對(duì)于公式(1)分別采用Crank-Nicholson格式(C-N格式)、上游加權(quán)法、人工擴(kuò)散量法，可得到以下三個(gè)差分方程：

(1) Crank-Nicholson格式的差分方程為

-ACi-1,n+1+BCi,n+1-DCi+1,n+1=ACi-1,n+

(4-B)Ci,n+DCi+1,n

(2)

式中：i為空間距離x的剖分離散結(jié)點(diǎn);n為時(shí)間t的剖分離散結(jié)點(diǎn);A、B、D為方程的系數(shù)，具體取值如下：

(2) 上游加權(quán)法的差分方程為

(3)

由上式可以發(fā)現(xiàn)，當(dāng)其他參數(shù)保持不變，u越大時(shí)ω取值要求越高。

(3) 人工擴(kuò)散量法的差分方程為

-ACi-1,n+1+BCi,n+1-DCi+1,n+1=ACi-1,n+

(4-B)Ci,n+DCi+1,n

(4)

式中：A、B、D為方程的系數(shù)，具體取值如下：

所謂過(guò)量就是由于數(shù)值解的振動(dòng)，出現(xiàn)了超過(guò)最大濃度和小于最小濃度的計(jì)算值，違背了基本的物理意義。數(shù)值彌散則是由截?cái)嗾`差引起的，這里我們采用Pe=Δx/α來(lái)描述對(duì)流和彌散的相對(duì)大小，Δx為空間步長(zhǎng)，α為縱向彌散率。Peclet數(shù)能夠有效地判定給定的網(wǎng)格剖分是否會(huì)引起濃度振蕩，對(duì)網(wǎng)格剖分具有指導(dǎo)意義[19]。按照Peclet數(shù)可以將對(duì)流-彌散方程分成兩大類：當(dāng)Peclet數(shù)小于或等于2時(shí)，為彌散占優(yōu)方程；當(dāng)Peclet大于2時(shí)，為對(duì)流占優(yōu)方程。

2三種差分格式的誤差及穩(wěn)定性分析

2.1Crank-Nicholson格式

采用泰勒級(jí)數(shù)展開(kāi)法分析可知，由于Crank-Nicholson格式空間項(xiàng)采用的是中心差商，所以其截?cái)嗾`差為O(Δt2)+O(Δx2)，下面來(lái)分析其穩(wěn)定性。

首先擴(kuò)充函數(shù)的定義域，使其在整個(gè)實(shí)軸上有定義，令

U(x,tn+1)=Ci,n+1

所以公式(2)可以改寫(xiě)為-AU(x-Δx,tn+1)+BU(x,tn+1)-DU(x+Δx,tn+1)=

AU(x-Δx,tn)+(4-B)U(x,tn)+DU(x+Δx,tn)

令

1-cos(kΔx)≥0,因?yàn)樯鲜降姆帜笧檎杂衸G|2-1≤0。

得出Crank-Nicholson格式是無(wú)條件穩(wěn)定的。

Crank-Nicholson格式是隱式的，即

F2(Ci+1,n+1,Ci,n+1,Ci-1,n+1)=F1(Ci+1,n,Ci,n,Ci-1,n)

Ci+1,n,Ci,n,Ci-1,n和兩端函數(shù)已知，只有Ci+1,n+1,Ci,n+1,Ci-1,n+1待求。因?yàn)楦咚?塞德?tīng)柕╗20]的收斂速度比較快，所以本文采用高斯-塞德?tīng)柕▉?lái)求解方程組。在此種迭代方法下，為了滿足迭代的收斂性，系數(shù)矩陣必須是嚴(yán)格對(duì)角占優(yōu)的，即|B|>|A|+|D|。將系數(shù)代入,并化簡(jiǎn)得

(5)

因此求解方程組時(shí)采用高斯-塞德?tīng)柗ǖ?，要求時(shí)間和空間剖分Δt、Δx以及彌散系數(shù)DL和實(shí)際速度u滿足式(5)Crank-Nicholson格式收斂。這里討論的計(jì)算時(shí)迭代的收斂性分析和差分格式的數(shù)值解對(duì)于精確解的收斂性是不同的概念。由于在建立差分方程以后，給定初值和邊界條件后，就必須通過(guò)解方程組來(lái)計(jì)算每個(gè)點(diǎn)的數(shù)值，本文在解方程時(shí)用高斯-塞德?tīng)柕?，所以在?jì)算時(shí)就有一個(gè)迭代過(guò)程的收斂性問(wèn)題，在文中“計(jì)算時(shí)迭代的收斂性”是指方程組的數(shù)值解是否收斂到方程組的精確解。另外兩種方法的收斂性也表示的是同樣的含義。

2.2上游加權(quán)法

由于公式(3)取的不是中心差分，所以截?cái)嗾`差是一階的，即為O(Δt)+O(Δx)。上游加權(quán)法的穩(wěn)定性分析與Crank-Nicholson格式類似，首先來(lái)擴(kuò)充函數(shù)的定義域，使其在整個(gè)實(shí)軸上有定義，每一項(xiàng)用Fourier積分來(lái)表示，引入歐拉公式，得出過(guò)度因子為G=m(4)/{-m(1)cos(kΔx)+m(2)-m(3)cos(kΔx)+i[m(1)sin(kΔx)-m(3)sin(kΔx)]}其中:

將系數(shù)代入上式，并整理得

H2[1+cos2(kΔx)]-2H2cos(kΔx)=

H2[1+cos2(kΔx)-2cos(kΔx)]≥0;

當(dāng)H≥0時(shí)，過(guò)度因子等式右端滿足小于或等于1，故增長(zhǎng)矩陣一致有界，式(3)是穩(wěn)定的，因此上游加權(quán)法的穩(wěn)定性條件是滿足下式：

(6)

由于差分格式是隱式的，與Crank-Nicholson格式類似，所以方程組采用高斯-塞德?tīng)柗ǖ鷷r(shí)，為了滿足迭代的收斂性，系數(shù)矩陣必須是嚴(yán)格對(duì)角占優(yōu)的,即

當(dāng)認(rèn)為u>0時(shí)，化簡(jiǎn)為

(7)

因此求解方程組時(shí)采用高斯-塞德?tīng)柗ǖ?，要求時(shí)間和空間剖分Δt、Δx以及彌散系數(shù)DL和實(shí)際速度u滿足式(7)時(shí)結(jié)果收斂。通過(guò)前面分析可知上游加權(quán)法的截?cái)嗾`差為O(Δt)+O(Δx)，說(shuō)明此方法與Crank-Nicholson格式相比解的精度降低了，尤其在速度很大時(shí)由(7)式收斂條件可知上游加權(quán)法就達(dá)不到收斂的要求。

2.3人工擴(kuò)散量法

(8)

用泰勒公式展開(kāi)上式的C(x+Δx)與C(x-Δx)，并整理得

說(shuō)明方程式(8)右端兩項(xiàng)都很準(zhǔn)確，再分析左端的時(shí)間項(xiàng)，對(duì)時(shí)間項(xiàng)進(jìn)行泰勒展開(kāi)并代入方程式(8)左端，有

所以式(1)真正離散后的方程變?yōu)?/p>

可見(jiàn)，為了校正這種誤差，引入人工擴(kuò)散系數(shù)，就此推出了人工擴(kuò)散量方法的表達(dá)式。

(9)

同理，若求解方程組時(shí)采用高斯-賽德?tīng)柗ǖ髸r(shí)間和空間剖分Δt、Δx以及彌散系數(shù)DL和〗實(shí)際速度u滿足公式(9)人工擴(kuò)散量法收斂。

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

本文考慮下述地下水溶質(zhì)運(yùn)移問(wèn)題，即半無(wú)限長(zhǎng)均質(zhì)砂柱，地下水以實(shí)際平均流速u穩(wěn)定流動(dòng)，初始濃度為0，其一端為定濃度邊界，假設(shè)無(wú)窮遠(yuǎn)處濃度分布為零，求其濃度分布。該問(wèn)題的解析解為

(10)

為了用無(wú)限長(zhǎng)的解析解近似有限長(zhǎng)的實(shí)際問(wèn)題，縮小誤差，本文通過(guò)多次試驗(yàn)確保數(shù)值解與解析解在有限長(zhǎng)的邊界處附近濃度都為0，取時(shí)間為4 d，空間為60 m，并取參數(shù)C0=1 mg/L、u=6 m/d、Δt=0.5 d、Δt=0.1 d、Δx=2 m，這里取DL=uα[21]，同時(shí)按照Pe=Δx/α取不同的彌散率來(lái)代表不同的對(duì)流強(qiáng)度，α取值分別為2 m、0.5 m、0.125 m和0.062 5 m。通過(guò)計(jì)算驗(yàn)證表明其都滿足三種差分格式的穩(wěn)定性和收斂性條件，對(duì)應(yīng)的對(duì)流占優(yōu)強(qiáng)度分別為1、4、16和32。圖1和圖2為對(duì)流強(qiáng)度分別為1和4時(shí)數(shù)值解與解析解的濃度對(duì)比圖，其中圖1和圖2中(a)和(b)分別為取時(shí)間段為1 d和3 d時(shí)三種差分格式隨空間位置的濃度分布圖。

圖1　當(dāng)對(duì)流強(qiáng)度為1時(shí)數(shù)值解與解析解的濃度對(duì)比圖Fig.1　Comparison of numerical and analytical solutions　　　for the concentration when Pe is 1

通過(guò)對(duì)比分析圖1和圖2可見(jiàn)，當(dāng)對(duì)流占優(yōu)較弱時(shí)，三種方法都很適用，其數(shù)值解與解析解幾乎重合，差分格式誤差都很??；當(dāng)對(duì)流強(qiáng)度為4時(shí)C-N格式開(kāi)始出現(xiàn)過(guò)量現(xiàn)象，上游加權(quán)法差分格式誤差最大。通過(guò)將圖2(b)與圖2(c)進(jìn)行對(duì)比分析可見(jiàn)，縮小時(shí)間步長(zhǎng)，三種方法的誤差都相對(duì)減小。

圖3和圖4為對(duì)流強(qiáng)度分別為16和32時(shí)數(shù)值解與解析解的濃度對(duì)比圖，其中圖3和圖4中(a)和(b)為取時(shí)間為3 d、時(shí)間步長(zhǎng)分別為0.5和0.1時(shí)三種差分格式隨空間位置的濃度分布圖。

通過(guò)對(duì)比分析圖3(a)和圖3(b)可見(jiàn)，人工擴(kuò)散量法對(duì)時(shí)間步長(zhǎng)比較敏感，時(shí)間步長(zhǎng)縮小時(shí)，差分格式誤差顯著變??；相比之下上游加權(quán)法在縮小時(shí)間步長(zhǎng)時(shí)，差分格式誤差變化并不明顯。

圖3　當(dāng)對(duì)流強(qiáng)度為16時(shí)數(shù)值解與解析解的濃度對(duì)比圖Fig.3　Comparison of numerical and analytical solutions　　　for the concentration when Pe is 16

圖4　當(dāng)對(duì)流強(qiáng)度為32時(shí)數(shù)值解與解析解的濃度對(duì)比圖Fig.4　Comparison of numerical and analytical solutions　　　for the concentration when Pe is 32

通過(guò)對(duì)比分析圖3和圖1及圖2可見(jiàn)，隨著對(duì)流強(qiáng)度的增大，C-N格式的過(guò)量現(xiàn)象越發(fā)明顯，數(shù)值解與解析解的總體誤差有所增大，這時(shí)人工擴(kuò)散量法數(shù)值解與解析解的誤差相比較小，尤其是在時(shí)間步長(zhǎng)較小時(shí)，此種方法開(kāi)始發(fā)揮優(yōu)勢(shì)，上游加權(quán)法精度不如人工擴(kuò)散量法。

通過(guò)對(duì)比分析圖4(a)和圖4(b)可見(jiàn)，當(dāng)對(duì)流強(qiáng)度為32時(shí)，人工擴(kuò)散量法在時(shí)間步長(zhǎng)較大時(shí)沒(méi)有出現(xiàn)過(guò)量現(xiàn)象，但在時(shí)間步長(zhǎng)過(guò)小時(shí)也出現(xiàn)了過(guò)量現(xiàn)象，說(shuō)明了在減小數(shù)值彌散的同時(shí)增大了解的振蕩；只有上游加權(quán)法沒(méi)有出現(xiàn)過(guò)量現(xiàn)象，但是誤差稍大，因?yàn)樯嫌渭訖?quán)法對(duì)空間步長(zhǎng)比較敏感，這時(shí)可以縮小空間步長(zhǎng)，通過(guò)更小的網(wǎng)格剖分，來(lái)減弱對(duì)流強(qiáng)度，以減小差分格式誤差。

4結(jié)論

采用數(shù)值法求解對(duì)流占優(yōu)的地下水溶質(zhì)運(yùn)移問(wèn)題時(shí)，難點(diǎn)是克服求解過(guò)程中的虛假數(shù)值現(xiàn)象。本文首先分析了差分格式的截?cái)嗾`差、穩(wěn)定性和收斂性，再進(jìn)行數(shù)值試驗(yàn)，結(jié)果表明:①在對(duì)流占優(yōu)比較弱時(shí)，三種方法都適用，但C-N格式最接近理論值；②當(dāng)對(duì)流強(qiáng)度增大時(shí)，例如數(shù)值試驗(yàn)中對(duì)流強(qiáng)度大于4時(shí)，建議使用時(shí)間步長(zhǎng)較小的人工擴(kuò)散量法，因?yàn)镃-N格式開(kāi)始出現(xiàn)過(guò)量現(xiàn)象，其誤差開(kāi)始增大，而上游加權(quán)法差分格式的誤差比前兩種方法稍大；③當(dāng)對(duì)流強(qiáng)度較大時(shí)，例如數(shù)值試驗(yàn)中對(duì)流強(qiáng)度大于32時(shí)，建議使用時(shí)間步長(zhǎng)較大的人工擴(kuò)散量法，或者使用空間步長(zhǎng)較小的上游加權(quán)法，通過(guò)調(diào)節(jié)空間步長(zhǎng)來(lái)減弱對(duì)流強(qiáng)度，再進(jìn)行求解。

對(duì)于一個(gè)實(shí)際問(wèn)題，彌散系數(shù)和流速是一定的，只有網(wǎng)格剖分可以調(diào)節(jié)，從前面分析可知當(dāng)彌散系數(shù)和流速已知時(shí)，可以適當(dāng)調(diào)節(jié)網(wǎng)格剖分，在滿足穩(wěn)定性和收斂性條件下，根據(jù)實(shí)際問(wèn)題在保證適當(dāng)計(jì)算量的前提下提高算法的精度。

參考文獻(xiàn)：

[1] 王焰新.地下水污染與防治[M].北京：高等教育出版社,2007.

[2] 盧曉華.基于數(shù)值模擬的企業(yè)地下水重金屬污染的環(huán)境影響預(yù)測(cè)評(píng)價(jià)[J].安全與環(huán)境工程,2014,21(1):93-97.

[3] 李合蓮,陳家軍.鈾尾礦庫(kù)中238U運(yùn)移數(shù)值模擬[J].安全與環(huán)境工程,2007,14(1):36-39.

[4] 鄧鼎興.某礦山地下水污染分區(qū)治理效果的數(shù)值模擬[J].安全與環(huán)境工程,2013,20(4):27-31.

[5] 魏恒,肖洪浪.地下水溶質(zhì)運(yùn)移模擬研究進(jìn)展[J].冰川凍土,2013,35(6):1584-1586.

[6] Hossain M A,Barber M E.Optimized Petrov-Galerkin model for advective-dispersive transport[J].AppliedMathematicsandComputation,2000,115(1):1-10.

[7] Dehghan M.Weighted finite difference techniques for the one-dimensional advection-diffusion equation[J].AppliedMathematicsandComputation,2004,147(2):307-319.

[8] Sun N Z,Yeh W G.A proposed upstream weight numerical method for simulation pollutant transporting ground water[J].WaterResourcesResearch,1983,19(6):1489-1500.

[9] Brooks A N,Hughes T J R.Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[J].ComputerMethodsinAppliedMechanicsandEngineering,1982,32(1/2/3):199-259.

[10]梅一,吳吉春.地下水溶質(zhì)運(yùn)移數(shù)值模擬中減少誤差的新方法[J].水科學(xué)進(jìn)展,2009,20(5):640-645.

[11]Douglas J,Russell T F.Numerical method for convection-dominated diffusion problems based on combing the method of characteristics with finite element or finite difference procedures[J].SIAMJournalonNumericalAnalysis,1982,19(5):871-885.

[12]Tsai T L,Yang J C,Huang L H.Characteristics method using cubic-spline interpolation for advection-diffusion equation[J].JournalofHydraulicEngineering,2004,130(6):580-585.

[13]Leismann H M,Frind E O.A symmetric-matrix time integration scheme for the efficient solution of advection-dispersion problems[J].WaterResourcesResearch,1989,25(6):1133-1139.

[14]劉揚(yáng).對(duì)流-擴(kuò)散方程的新型Crank-Nicholson差分格式[J].數(shù)學(xué)雜志,2005,25(4):464-467.

[15]成建梅,胡進(jìn)武.飽和水流溶質(zhì)運(yùn)移問(wèn)題數(shù)值解法綜述[J].水文地質(zhì)工程地質(zhì),2003,30(2):99-101.

[16]劉苑,武曉峰.地下水中污染物運(yùn)移過(guò)程數(shù)值模擬算法的比較[J].環(huán)境工程學(xué)報(bào),2008,2(2):230-234.

[17]寧一鳴.地下水溶質(zhì)運(yùn)移對(duì)流-彌散方程數(shù)值解方法對(duì)比研究[D].武漢:中國(guó)地質(zhì)大學(xué),2008.

[18]陳崇希,李國(guó)敏.地下水溶質(zhì)運(yùn)移理論及模型[M].武漢:中國(guó)地質(zhì)大學(xué)出版社,1996.

[19]馬倩,李海龍.變異Henry問(wèn)題檢驗(yàn)網(wǎng)格Peclet數(shù)在數(shù)值計(jì)算中的重要性[J].水文地質(zhì)工程地質(zhì),2014,41(3):1-5.

[20]李慶揚(yáng),王能超,易大義.數(shù)值分析[M].北京:清華大學(xué)出版社,2008.

[21]Fetter C W.ContaminantHydrology[M].New Jersey:Prentice-Hal,1999:45-186.

Empirical Study of the Finite Difference Schemes of Solute Transport Equation in Groundwater

JIAO Tian1,WANG Junxia2,TANG Zhonghua1,LIU Geng3

(1.SchoolofEnvironmentalStudies,ChinaUniversityofGeosciences,Wuhan430074,China;2.SchoolofMathematicsandPhysics,ChinaUniversityofGeosciences,Wuhan430074,China;3.BeijingDeliheEnvironmentalGovernanceLtd.,Beijing100025,China)

Abstract:In solving the advection-dominated solute transport equation in groundwater,in order to overcome the excess and numerical dispersion produced by finite-difference and finite element methods,many numerical methods are adopted,such as Crank-Nicholson scheme,upstream weighting method and artificial diffusion method.Taking the one-dimensional groundwater solute transport equation as an example,this paper analyses the truncation error,stability and convergence of the three difference schemes,and compares the numerical solutions of different algorithms with the analytical solutions in different advection-dominated intensities through numerical experiments in order to explore the advantages and disadvantages and applicable conditions of different algorithms.The results indicate that when the advection-dominated intensity is weak,all the methods are applicable,while the error of CN scheme is the minimum;with the increasing advection-dominated intensity,the artificial diffusion method with smaller time step is applicable;once the advection-dominated intensity is relatively strong,the artificial diffusion method with larger time step or the upstream weighting method with smaller space step is applicable.This research provides a basis for the selection of algorithm for solving the practical problems.

Key words:groundwater solute transport equation;advection-dominated intensity;finite difference scheme;analytical solution;empirical comparison

文章編號(hào)：1671-1556(2016)03-0017-07

收稿日期：2015-10-11修回日期：2016-03-20

基金項(xiàng)目：中國(guó)博士后科學(xué)基金項(xiàng)目(201404406064);中國(guó)地質(zhì)調(diào)查局項(xiàng)目(1212011121142)

作者簡(jiǎn)介：焦甜(1993—)，女，碩士研究生，主要研究方向?yàn)榈叵滤叭苜|(zhì)運(yùn)移數(shù)值模擬。E-mail:jiaotianedu@sina.cn

中圖分類號(hào)：X143;P641.2

文獻(xiàn)標(biāo)識(shí)碼：A

DOI：10.13578/j.cnki.issn.1671-1556.2016.03.003

通訊作者：唐仲華(1958—)，男，教授，主要從事滲流理論及其數(shù)值模擬方面的研究。E-mail:zhhtang@cug.edu.cn