蔡松柏，邵明安，齊麗彬

(1.湖南師范大學(xué) 資源與環(huán)境科學(xué)學(xué)院,湖南 長(zhǎng)沙 410081;2.中國(guó)科學(xué)院 水利部水土保持研究所,陜西 楊凌 712100;3.三門(mén)峽職業(yè)技術(shù)學(xué)院 生化工程系,河南 三門(mén)峽 472000)

由于工農(nóng)業(yè)生產(chǎn)、雨水徑流沖刷和人為因素等作用，造成水流中溶質(zhì)瞬間注入或連續(xù)注入從而引起河流或地下水的水質(zhì)改變。研究溶質(zhì)瞬時(shí)注入或連續(xù)注入情況下不同時(shí)刻溶質(zhì)濃度在河流或地下水流中的分布狀況、及不同位置時(shí)溶質(zhì)濃度隨時(shí)間的變化情況，可通過(guò)數(shù)學(xué)模型來(lái)模擬或預(yù)報(bào)水質(zhì)在時(shí)間或空間的變化，達(dá)到為水環(huán)境質(zhì)量預(yù)報(bào)、控制和水資源管理服務(wù)的目的[1-22]。

數(shù)值模擬水流溶質(zhì)運(yùn)移的主要方法有有限元法、邊界單元法、有限差分法和其他數(shù)值法等。由于進(jìn)行水質(zhì)模擬過(guò)程中涉及的參數(shù)在不同的時(shí)間和空間上變化較大，在水質(zhì)模擬過(guò)程中，會(huì)造成數(shù)值彌散、數(shù)值波動(dòng)和計(jì)算量巨大等問(wèn)題。國(guó)內(nèi)主要的研究有徐玉佩[2-3]用Galerkin法推導(dǎo)了二維地下水運(yùn)動(dòng)的有限元方程和溶質(zhì)運(yùn)移有限元方程，并編制了計(jì)算程序應(yīng)用于明溝排水和田面保持淹灌水層情況下的沖洗種稻改良鹽堿地的分析,還與現(xiàn)有解析解和實(shí)驗(yàn)資料進(jìn)行了對(duì)照。徐玉佩還用水動(dòng)力彌散理論進(jìn)行了地下水中水鹽動(dòng)態(tài)的研究,可直接用于鹽堿地的預(yù)防和改良。黃康樂(lè)[4]則提出了一種求解二維飽和-非飽和溶質(zhì)運(yùn)移問(wèn)題的交替方向特征有限單元法(ADCG)，在克服數(shù)值彌散、數(shù)值波動(dòng)和提高計(jì)算速度方面具有顯著的優(yōu)越性。張效先等[5]提出了一種計(jì)算地下水中溶質(zhì)運(yùn)移的數(shù)值法，將多維問(wèn)題分解成幾個(gè)一維問(wèn)題, 減少了計(jì)算工作量, 并消除了數(shù)值波動(dòng)現(xiàn)象,具有較高的精度。馮紹元[6]提出了特征有限元法的基本理論和實(shí)施過(guò)程,其計(jì)算實(shí)例表明能消除大部分的數(shù)值彌散而又不增加解的振蕩性。姚磊華[7]首先基于分步求解的思想, 用廣義迎風(fēng)對(duì)偶單元均衡法求解對(duì)流定解問(wèn)題, 對(duì)擴(kuò)散定解問(wèn)題采用一般的Galerkin有限元法求解, 不僅避免了數(shù)值彌散和過(guò)量問(wèn)題, 而且避免了求節(jié)點(diǎn)速度這一步, 簡(jiǎn)化了運(yùn)算步驟。邱克儉等[8]結(jié)合拉格朗日觀點(diǎn)的運(yùn)動(dòng)坐標(biāo)與歐拉觀點(diǎn)的固定坐標(biāo),分別求解對(duì)流-彌散方程的對(duì)流效應(yīng)和彌散效應(yīng), 并通過(guò)與理論解和試驗(yàn)結(jié)果做比較，驗(yàn)證了模型的可靠性。任理等[9]采用混合拉普拉斯變換有限單元法求解對(duì)流占優(yōu)的地下水溶質(zhì)運(yùn)移問(wèn)題，能有效地消除數(shù)值擴(kuò)散和過(guò)量現(xiàn)象，特別適于大區(qū)域地下水污染的長(zhǎng)期預(yù)報(bào)。李煥榮等[10-11]利用廣義差分法建立了一維非飽和水流問(wèn)題的守恒形式的數(shù)值模型, 具有計(jì)算量小和穩(wěn)定性好的特點(diǎn), 其提出的非粘性土壤水中溶質(zhì)運(yùn)移問(wèn)題的守恒混合元格式，精度高且數(shù)值穩(wěn)定。段德宏等[14]以黃河蘭州段為例,采用一維穩(wěn)態(tài)河流水質(zhì)模型,對(duì)岸邊19個(gè)污水排放口及共約32.9 km長(zhǎng)的全河段污染狀況成功地進(jìn)行了模擬。徐文彬等[15]根據(jù)數(shù)值隨時(shí)間步長(zhǎng)變化的特點(diǎn)，提出了針對(duì)數(shù)值波動(dòng)、彌散和過(guò)量現(xiàn)象的新改進(jìn)方案，并用實(shí)例驗(yàn)證了方法的可行性。何丕文[16]利用一維水質(zhì)模型，采用有限差分法求解河流中污染物濃度，以淮河淮南段為例進(jìn)行了驗(yàn)證。馬東豪等[17]研究了土壤溶質(zhì)遷移的兩流區(qū)與兩區(qū)模型，結(jié)果表明與兩區(qū)模型相比,兩流區(qū)模型可以更好地描述優(yōu)先流情況下的溶質(zhì)穿透曲線。對(duì)于無(wú)優(yōu)先流情況, 連續(xù)流輸入情況下兩個(gè)模型的參數(shù)一致性較好；而在脈沖輸入時(shí)一致性較差。國(guó)外在溶質(zhì)運(yùn)移模型及其數(shù)值計(jì)算方法上的研究也很多，如JAISWAL等[18]運(yùn)用Chebyshev配點(diǎn)法求解一維溶質(zhì)運(yùn)移問(wèn)題的數(shù)值解。SILAVWE等[19]討論了某些數(shù)值方法用于一維溶質(zhì)運(yùn)移模型的參數(shù)估計(jì)，通過(guò)將數(shù)值解與試驗(yàn)結(jié)果比較，反演溶質(zhì)運(yùn)移模型參數(shù)。JAISWAL等[20]采用算子矩陣法對(duì)多孔介質(zhì)的非線性微分方程進(jìn)行了數(shù)值分析，該法應(yīng)用切比雪夫多項(xiàng)式并結(jié)合切比雪夫?qū)?shù)運(yùn)算矩陣和譜配置法對(duì)解進(jìn)行逼近。其優(yōu)點(diǎn)是可以將這些問(wèn)題轉(zhuǎn)化為易于求解的代數(shù)方程組。CARR[21]利用拉普拉斯變換導(dǎo)出了層狀多孔介質(zhì)溶質(zhì)運(yùn)移的新半解析解，該法在相鄰層的界面引入表示運(yùn)移量的未知函數(shù)，將多層問(wèn)題置于拉普拉斯域上分層解決，然后再將其數(shù)值反演回時(shí)域。ALLWRIGHT等[22]對(duì)地下水溶質(zhì)運(yùn)移方程在對(duì)流明顯強(qiáng)于擴(kuò)散時(shí)出現(xiàn)的數(shù)值波動(dòng)提出了兩種新的數(shù)值近似計(jì)算方法，即只適用于對(duì)流項(xiàng)的迎風(fēng)Crank-Nicolson格式和加權(quán)迎風(fēng)-下風(fēng)格式，并對(duì)這些新提出的格式進(jìn)行了數(shù)值穩(wěn)定性分析和比較。發(fā)現(xiàn)如果Crank-Nicolson格式只用于對(duì)流項(xiàng)時(shí)很合宜。此外，顯式加權(quán)迎風(fēng)-下風(fēng)有限差分格式是對(duì)傳統(tǒng)顯式一階迎風(fēng)格式的改進(jìn)，而隱式加權(quán)一階迎風(fēng)-下風(fēng)有限差分格式在給定適當(dāng)?shù)臋?quán)重因子時(shí)無(wú)條件穩(wěn)定。

本文主要貢獻(xiàn)在于對(duì)一維溶質(zhì)對(duì)流彌散基本方程，首先，建立了一維溶質(zhì)運(yùn)移3種變?nèi)苜|(zhì)運(yùn)移參數(shù)模型，即線性變化模型/多項(xiàng)式變化模型、指數(shù)變化模型和冪變化模型；其次，采用顯式時(shí)域差分法成功地進(jìn)行了溶質(zhì)濃度時(shí)空分布的數(shù)值計(jì)算，有效克服了很多研究人員遭遇的數(shù)值波動(dòng)、數(shù)值彌散和計(jì)算量大等難題，并與解析解和實(shí)驗(yàn)結(jié)果做了對(duì)比，驗(yàn)證了該數(shù)值方法的正確性；最后，對(duì)3種溶質(zhì)運(yùn)移參數(shù)模型進(jìn)行了參數(shù)分析，獲得了3種模型的初步認(rèn)識(shí)。

1 基本原理

研究溶質(zhì)在水流中的遷移規(guī)律，試驗(yàn)測(cè)定可采用電解質(zhì)脈沖示蹤法等來(lái)進(jìn)行。電解質(zhì)脈沖在水流中的運(yùn)移規(guī)律符合溶質(zhì)對(duì)流彌散基本方程，根據(jù)夏衛(wèi)生等[12-13]的研究，其溶質(zhì)運(yùn)移參數(shù)隨水流位置而變化，在大規(guī)模的地下水流問(wèn)題中，溶質(zhì)運(yùn)移參數(shù)亦會(huì)隨空間位置而變化。對(duì)其溶質(zhì)運(yùn)移規(guī)律的模擬可采用近似解析解和數(shù)值解。

1.1 一維變?nèi)苜|(zhì)運(yùn)移參數(shù)對(duì)流彌散溶質(zhì)運(yùn)移基本方程及定解問(wèn)題

由邵明安等[1]提出的一維變?nèi)苜|(zhì)運(yùn)移參數(shù)對(duì)流彌散溶質(zhì)運(yùn)移基本方程可寫(xiě)為：

(1)

在脈沖注入時(shí)的初始和邊界條件可寫(xiě)為：

c(x,0)=c0

(2)

c(0,t≤t0)=c1

(3)

c(0,t>t0)=c0

(4)

c(l,t)=c0

(5)

式中：x為水流方向長(zhǎng)度坐標(biāo)，t為時(shí)間，c(x,t)為水流中溶質(zhì)濃度，l為計(jì)算水流長(zhǎng)度，t0為溶質(zhì)脈沖注入的時(shí)間，c1為脈沖注入溶質(zhì)的濃度，c0為水流中溶質(zhì)的初始和平衡濃度，D(x)為水流中溶質(zhì)擴(kuò)散系數(shù)，u(x)為水流速度。本文假定D(x)，u(x)為下列3種模式。

1)線性模式/多項(xiàng)式模式

即,線性模式：

D(x)=D0(1+bx/l)

(6)

u(x)=u0(1+bux/l)

(7)

和多項(xiàng)式模式：

D(x)=D0(1+bx/l+cx2/l2)

(8)

u(x)=u0(1+bux/l+cux2/l2)

(9)

2)冪函數(shù)模式

D(x)=D0(1+cxb)

(10)

u(x)=u0(1+cuxbu)

(11)

3)指數(shù)函數(shù)模式

D(x)=D0(1+cbx)

(12)

(13)

式中：b，c，bu和cu為模型參數(shù)。

1.2 一維常溶質(zhì)運(yùn)移參數(shù)溶質(zhì)對(duì)流彌散基本方程的解析解

定義溶質(zhì)擴(kuò)散系數(shù)和水流速度的平均值分別為：

(14)

(15)

對(duì)于式(6)和(7)線性模型，則有：

(16)

(17)

由文獻(xiàn)[1]，利用Heavyside階躍函數(shù)，本文首次得到了D(x)=D，u(x)=u時(shí)式(1)的解析解為：

(18)

2 理論與列式

2.1 差分網(wǎng)格

由于溶質(zhì)濃度隨時(shí)間和空間連續(xù)變化，可沿x方向取水流有限長(zhǎng)度l作為計(jì)算對(duì)象，將l分為等間距的n個(gè)單元，結(jié)點(diǎn)編號(hào)為i，i=0,1,2,…,n，步長(zhǎng)為Δx=l/n，將時(shí)間坐標(biāo)劃分為時(shí)間步長(zhǎng)為Δt的時(shí)段，結(jié)點(diǎn)編號(hào)為k，k=0,1,2,…,m。

2.2 即時(shí)顯式差分格式

(19)

將式(19)變形，可得差分遞進(jìn)格式：

i=1,2,3,…,n-1

(20)

其中：

(21)

(22)

(23)

3 計(jì)算結(jié)果與討論

根據(jù)上面的理論基礎(chǔ)，利用雙精度FORTRAN語(yǔ)言編寫(xiě)程序，分別采用方程(19)所代表的數(shù)值方法和方程(18)所代表的解析方法求得溶質(zhì)濃度分布，以檢驗(yàn)數(shù)值方法的可靠性。并對(duì)不同溶質(zhì)運(yùn)移參數(shù)變化模式下瞬時(shí)注入溶質(zhì)后水流中溶質(zhì)濃度的時(shí)空分布進(jìn)行了數(shù)值模擬和參數(shù)分析。

3.1 一維溶質(zhì)運(yùn)移參數(shù)線性變化時(shí)方程數(shù)值解與解析解比較

由夏衛(wèi)生等[12-13]的試驗(yàn)結(jié)果得知，坡面薄層水流的溶質(zhì)擴(kuò)散系數(shù)及水流速度均隨位置而變化。本計(jì)算的變?nèi)苜|(zhì)運(yùn)移參數(shù)模型采用D(x)=D0(1+bx/l)，u(x)=u0(1+bux/l)。根據(jù)夏衛(wèi)生等[12-13]的試驗(yàn)數(shù)據(jù)，取D0=5.1×10-2m2/s，b=-0.55，u0=1.02 m/s和bu=0.042。l=3.5 m，t0=0.6 s， 注入溶液濃度c1=0.63，水流起始溶質(zhì)濃度和平衡濃度c0=0.20，Δt=0.04 sec和Δx=0.025 m，分別求得了水流方向上不同位置處溶質(zhì)濃度的數(shù)值解和近似解析解。

圖1給出了數(shù)值解、近似解析解與文獻(xiàn)[12]實(shí)測(cè)結(jié)果的比較，從圖中可以看出兩種模擬結(jié)果同夏衛(wèi)生等[12-13]的結(jié)果一致，測(cè)量曲線和模擬曲線的溶質(zhì)濃度最大值均呈逐漸減小的趨勢(shì)；在距離較大處其模擬結(jié)果和實(shí)測(cè)結(jié)果吻合較好。近似解析解和數(shù)值解能較好地吻合，特別是距離越大，二者吻合得越好。這說(shuō)明顯式差分法求解溶質(zhì)運(yùn)移方程的精度已經(jīng)能夠達(dá)到要求。

圖1 模擬結(jié)果和實(shí)測(cè)結(jié)果比較

表1給出了一維常溶質(zhì)運(yùn)移參數(shù)模型在脈沖注入后溶質(zhì)濃度隨深度、時(shí)間變化的數(shù)值解式(20)與解析解式(18)的比較，其中，常溶質(zhì)運(yùn)移參數(shù)D0=5.1×10-2m2/s，u0=1.02 m/s，其余計(jì)算參數(shù)同圖1取值，從表中可以看出二者吻合非常好，說(shuō)明了數(shù)值方法是可靠的。

表1 常系數(shù)模型注入脈沖溶質(zhì)后溶質(zhì)濃度變化的數(shù)值解與解析解的比較

由圖1和表1可知，常溶質(zhì)運(yùn)移參數(shù)一維溶質(zhì)運(yùn)移基本方程的解析解與變?nèi)苜|(zhì)運(yùn)移參數(shù)一維溶質(zhì)運(yùn)移基本方程的數(shù)值解基本一致。因此，求解溶質(zhì)運(yùn)移基本方程可以采用本文數(shù)值方法。

3.2 一維溶質(zhì)運(yùn)移參數(shù)多項(xiàng)式變化時(shí)的數(shù)值結(jié)果

給定D0=0.03 m2/s，u0=1.00 m/s，計(jì)算水流位置l=3.5 m，脈沖注入溶質(zhì)濃度時(shí)間t0=0.6 s，注入溶液濃度c1=1.00，水流起始溶質(zhì)濃度和平衡濃度c0=0.20，計(jì)算步長(zhǎng)Δt=0.04 sec和Δx=0.025 m,當(dāng)式(8)和(9)中模型參數(shù)c=0.1，b=±0.6,±0.4,±0.2，對(duì)cu,bu和c,b同步一樣取完全相同值時(shí)，其數(shù)值解法求出的溶質(zhì)濃度沿水流位置的分布規(guī)律如圖2所示，圖2中c和b同時(shí)代表cu,bu和c,b的取值。從圖2可以看出，c,cu的零值附近微小變化對(duì)溶質(zhì)濃度分布影響較計(jì)算方法帶來(lái)的差別要小。c,cu值一定時(shí)，隨著b,bu從-0.6到0.6的增加，某時(shí)刻溶質(zhì)濃度位置上的分布，其峰值和峰寬增加，波峰位置滯后。時(shí)間越長(zhǎng)，這種濃度分布的波峰變化規(guī)律越明顯。在1 s時(shí)間以?xún)?nèi)，本文對(duì)b,bu和c,cu的不同取值對(duì)溶質(zhì)濃度分布影響不大。

圖2 多項(xiàng)式模式溶質(zhì)參數(shù)的一維溶質(zhì)運(yùn)移濃度分布及比較

3.3 一維溶質(zhì)運(yùn)移參數(shù)冪函數(shù)變化時(shí)的數(shù)值結(jié)果

溶質(zhì)運(yùn)移參數(shù)采用式(10)和(11)的冪函數(shù)變化模型，給定模型參數(shù)b=2且c=±0.9,±0.6,±0.3下，cu,bu和c,b取值完全相同，其余計(jì)算參數(shù)D0，u0，l，t0，c1，c0，Δt和Δx同以上多項(xiàng)式模式中參數(shù)的取值，一維溶質(zhì)運(yùn)移濃度分布數(shù)值解如圖3所示。圖3表明，參數(shù)中指數(shù)b/c，bu/cu較大時(shí)，c,cu值的大小對(duì)溶質(zhì)濃度分布的影響較小，時(shí)間越短，這種影響越小。

圖3 冪函數(shù)模式溶質(zhì)參數(shù)的一維溶質(zhì)運(yùn)移濃度分布及比較

3.4 一維溶質(zhì)運(yùn)移參數(shù)指數(shù)變化時(shí)的數(shù)值結(jié)果

圖4所示為溶質(zhì)運(yùn)移參數(shù)模型采用式(12)和(13)模型的數(shù)值計(jì)算結(jié)果，計(jì)算中取模型參數(shù)b=0.125且c=±0.9,±0.6,±0.3下，cu,bu和c,b同步取完全相同的值，其余計(jì)算參數(shù)D0，u0，l，t0，c1，c0，Δt和Δx同以上多項(xiàng)式模式中參數(shù)的取值。從圖4看出模型參數(shù)的變化對(duì)濃度分布形態(tài)的影響較小，對(duì)于本例，當(dāng)b,bu一定時(shí)，c,cu的增大會(huì)導(dǎo)致濃度分布峰值位置滯后，濃度峰值增加。c,cu一定時(shí)，b,bu在一定范圍內(nèi)的變化對(duì)溶質(zhì)濃度分布的影響很小。

圖4 指數(shù)模式溶質(zhì)參數(shù)的一維溶質(zhì)運(yùn)移濃度分布及比較

4 結(jié)論

本文采用顯式時(shí)間差分格式計(jì)算變?nèi)苜|(zhì)運(yùn)移參數(shù)模型一維水流溶質(zhì)運(yùn)移過(guò)程，計(jì)算結(jié)果通過(guò)解析解進(jìn)行驗(yàn)證，并和試驗(yàn)結(jié)果進(jìn)行對(duì)比分析，說(shuō)明算法可靠，顯式法列式簡(jiǎn)單，更適于復(fù)雜模型計(jì)算[1]；模型越復(fù)雜，本法的計(jì)算效率越高。本文首次進(jìn)行了變?nèi)苜|(zhì)運(yùn)移參數(shù)模型溶質(zhì)運(yùn)移問(wèn)題的數(shù)值計(jì)算，并建立不同的溶質(zhì)運(yùn)移參數(shù)模型，對(duì)該模型進(jìn)行了參數(shù)分析。提出的系列算法和處理技巧可以應(yīng)用到溶質(zhì)注入水流后，變?nèi)苜|(zhì)運(yùn)移參數(shù)模型的溶質(zhì)遷移規(guī)律研究，通過(guò)建立不同模式下溶質(zhì)運(yùn)移參數(shù)模型，從而實(shí)現(xiàn)復(fù)雜條件下溶質(zhì)運(yùn)移過(guò)程的預(yù)測(cè)和模擬。