一維變參數(shù)溶質(zhì)運(yùn)移模型的數(shù)值研究

蔡松柏，邵明安，齊麗彬

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

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

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

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

1 基本原理

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

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

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

(1)

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

c(x,0)=c0

(2)

c(0,t≤t0)=c1

(3)

c(0,t>t0)=c0

(4)

c(l,t)=c0

(5)

式中：x為水流方向長度坐標(biāo)，t為時(shí)間，c(x,t)為水流中溶質(zhì)濃度，l為計(jì)算水流長度，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ì)對流彌散基本方程的解析解

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

(14)

(15)

對于式(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方向取水流有限長度l作為計(jì)算對象，將l分為等間距的n個(gè)單元，結(jié)點(diǎn)編號為i，i=0,1,2,…,n，步長為Δx=l/n，將時(shí)間坐標(biāo)劃分為時(shí)間步長為Δt的時(shí)段，結(jié)點(diǎn)編號為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語言編寫程序，分別采用方程(19)所代表的數(shù)值方法和方程(18)所代表的解析方法求得溶質(zhì)濃度分布，以檢驗(yàn)數(shù)值方法的可靠性。并對不同溶質(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í)測結(jié)果的比較，從圖中可以看出兩種模擬結(jié)果同夏衛(wèi)生等[12-13]的結(jié)果一致，測量曲線和模擬曲線的溶質(zhì)濃度最大值均呈逐漸減小的趨勢；在距離較大處其模擬結(jié)果和實(shí)測結(jié)果吻合較好。近似解析解和數(shù)值解能較好地吻合，特別是距離越大，二者吻合得越好。這說明顯式差分法求解溶質(zhì)運(yùn)移方程的精度已經(jīng)能夠達(dá)到要求。

圖1 模擬結(jié)果和實(shí)測結(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取值，從表中可以看出二者吻合非常好，說明了數(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ì)算步長Δt=0.04 sec和Δx=0.025 m,當(dāng)式(8)和(9)中模型參數(shù)c=0.1，b=±0.6,±0.4,±0.2，對cu,bu和c,b同步一樣取完全相同值時(shí)，其數(shù)值解法求出的溶質(zhì)濃度沿水流位置的分布規(guī)律如圖2所示，圖2中c和b同時(shí)代表cu,bu和c,b的取值。從圖2可以看出，c,cu的零值附近微小變化對溶質(zhì)濃度分布影響較計(jì)算方法帶來的差別要小。c,cu值一定時(shí)，隨著b,bu從-0.6到0.6的增加，某時(shí)刻溶質(zhì)濃度位置上的分布，其峰值和峰寬增加，波峰位置滯后。時(shí)間越長，這種濃度分布的波峰變化規(guī)律越明顯。在1 s時(shí)間以內(nèi)，本文對b,bu和c,cu的不同取值對溶質(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值的大小對溶質(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ù)的變化對濃度分布形態(tài)的影響較小，對于本例，當(dāng)b,bu一定時(shí)，c,cu的增大會導(dǎo)致濃度分布峰值位置滯后，濃度峰值增加。c,cu一定時(shí)，b,bu在一定范圍內(nèi)的變化對溶質(zhì)濃度分布的影響很小。

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

4 結(jié)論

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