一維對流彌散方程的顯式差分法求解及其收斂性分析

(石家莊鐵道大學(xué) 土木工程學(xué)院，河北 石家莊 050043)

1 研究背景

在濱海及海島地區(qū)，由于地下水的超量開采，海水入侵已經(jīng)成為沿海地區(qū)普遍面臨的重大地質(zhì)災(zāi)害之一。自20世紀(jì)70年代以來，國內(nèi)外對海水入侵進(jìn)行了大量研究并取得很大進(jìn)展，其中數(shù)學(xué)模型的研究是其中的重要方面之一，而對流彌散模型能夠較好地模擬和解釋海水入侵問題[1]。研究一維對流彌散方程是研究海水入侵的重要方法，能夠粗略地模擬出海水入侵情況，相比較二維和三維所需要的參數(shù)較少，是一種簡易可行的方法。楊金忠[2]對各種求解對流彌散方程數(shù)值方法比較，認(rèn)為對于一維問題，有限差分法能夠得到精確結(jié)果，因此本文提出利用顯式有限差分法對一維對流彌散模型進(jìn)行求解，具有操作簡單、穩(wěn)定性和收斂性好的特點，是求解一維對流彌散模型的有效方法。

有限差分法是一種古典的數(shù)值計算方法。隨著數(shù)值計算方法的研究和發(fā)展，差分法已經(jīng)被廣泛應(yīng)用在地下水滲流問題和濃度擴(kuò)散問題中。其基本思想是：用滲流區(qū)內(nèi)有限個離散點的集合代替連續(xù)的滲流區(qū)，在這些離散點上用差商近似地代替微商，將微分方程及其定解條件化為以未知函數(shù)在離散點上的近似值為未知量的代數(shù)方程(稱為差分方程)，然后求解差分方程，從而得到微分方程的解在離散點上的近似值。有限差分法難免會產(chǎn)生截斷誤差、舍入誤差和測量誤差，如果在某個節(jié)點某個時階處出現(xiàn)誤差后，會對下一步的迭代求解過程產(chǎn)生影響，是誤差逐漸變小還是被無限放大，需要討論差分法的收斂性和穩(wěn)定性[3-4]。

夏源等[5]提出分?jǐn)?shù)階對流彌散方程的數(shù)值求解方法，楊淑伶[6]提出求解分?jǐn)?shù)階對流彌散方程的邊界值法。傳統(tǒng)研究對流彌散方程的數(shù)值求解時，往往將縱向彌散系數(shù)看做常數(shù)，而事實上縱向彌散系數(shù)與流速的大小成正比，當(dāng)流速在空間的分布不相等時，流速和縱向彌散系數(shù)均不能看做常數(shù)。在研究對流彌散方程的差分法求解時，在不同滲流方向下，流速會與鹽分遷移的方向相同或相反，代入到差分式當(dāng)中會有正負(fù)號的影響，而彌散系數(shù)只與流速的大小有關(guān)，為保證迭代式中的流速項系數(shù)為正值，在不同的滲流方向下，對流彌散方程會出現(xiàn)2種差分迭代式，這是以往研究沒有考慮到的。

當(dāng)利用有限差分法求解一維對流彌散方程時，需要考慮邊界條件的問題，李新潔等[7]采用Dirichlet零邊值邊界條件求解對流彌散方程。本文引入邊界衰減因子來解釋鹽分在邊界處存在稀釋和衰減的情況，并討論在不同滲流方向下的影響，最后討論了在不同滲流方向下給出不同的差分式，并給出了各自對應(yīng)的收斂條件，在收斂條件下，利用顯式差分法求解對流彌散方程是收斂的。

2 研究方法

2.1 一維對流彌散方程

以一維對流彌散方程為研究對象，不計源匯項，且孔隙率n為常數(shù)，在穩(wěn)定滲流場中，若求解對流彌散方程得出溶質(zhì)濃度變化引起的水密度變化可以忽略，則水流方程和溶質(zhì)遷移方程可以獨立求解。由水流方程的解得出研究區(qū)域及時段的速度分量，然后把速度作為輸入代入對流彌散方程。這種“去耦”法計算效率高[8]。

(1)

(2)

Dl=αL|Vx|+D*

(3)

式中，k為滲透系數(shù)；Dl為彌散系數(shù)，為反映可溶性物質(zhì)通過滲透介質(zhì)時的彌散現(xiàn)象強(qiáng)弱的物理量；αL為縱向彌散度；D*=τD0，τ為孔隙介質(zhì)的彎曲度，取值范圍為0.5～0.01，D0為開放水體的分子擴(kuò)散系數(shù)，對于NaCl鹽溶液經(jīng)查表得D0=1.48×10-5cm2/s[9]，τ取0.5，代入得到D*=7.4×10-10m2/s，其值太小可忽略不計。

Dl=αL|Vx|

(4)

2.2 劃分網(wǎng)格

在海水入侵模型中，假設(shè)左側(cè)海水高度為h0，右側(cè)淡水水頭為hN，滲流區(qū)間[0，L]均勻分布N+1個節(jié)點，則空間步長Δ=L/N，X0=0，節(jié)點號數(shù)x乘以空間步長Δx便是該節(jié)點距左側(cè)海水邊界的長度，lxi=iΔ(=1,2，…，N)；取時間步長為Δt，則tn=nΔt(n=1,2, …,M)。

2.3 微分方程的顯式差分化

對于一維海水入侵對流彌散模型，存在2種類型，正向滲流類型是左側(cè)海水水位高，右側(cè)淡水水位低，即h0>hN,如圖1所示，此時速度方向與鹽分遷移的方向相同，且與x軸方向相同；逆向滲流類型是左側(cè)海水水位低，右側(cè)淡水水位高，即h0

在正向滲流的類型下，流速為正值，根據(jù)達(dá)西定律，待到穩(wěn)定后，Vi的值和hi存在一一對應(yīng)關(guān)系，根據(jù)空間節(jié)點處流量相等的原理[10]，可以求得穩(wěn)定后的流速表達(dá)式為

(5)

(6)

(7)

(8)

將一維對流彌散方程式(1)差分化為

(9)

(10)

對式(10)進(jìn)行恒等變換有

(11)

令

(12)

則可化簡為

(13)

式(13)代表一維對流彌散方程在正向滲流類型下的顯式差分迭代公式。

在逆向滲流的情況下，代入流速的計算值為負(fù)值，考慮到彌散系數(shù)不隨流速正負(fù)的影響而只與大小有關(guān)，為了保證代入的顯式差分式中的流速項系數(shù)為正值，故將式(1)差分化為

(14)

當(dāng)i=N時，濃度的差分式與式(10)相同，且此時的速度表達(dá)式為

(15)

(16)

同理化為標(biāo)準(zhǔn)的差分式為

(17)

式(17)代表一維對流彌散方程在逆向滲流類型下的顯式差分迭代公式。

2.4 差分格式的收斂性判定

前面建立差分方程時，是用差商代替微商，在數(shù)值法中提出一種差分格式，需要了解其收斂性。

記

(18)

2.4.1 正向滲流類型的收斂條件

當(dāng)為正向滲流類型時，左側(cè)海水水位高于右側(cè)淡水水位，根據(jù)式(10)得到

(19)

令

(20)

有

(21)

當(dāng)滿足

(22)

根據(jù)絕對值不等式有

(23)

當(dāng)?shù)谝环N類型下，Vi>Vi-1，故有

(24)

如此類推，對于某個時刻tM有

(25)

(26)

2.4.2 逆向滲流類型的收斂條件

當(dāng)為逆向滲流類型時，左側(cè)海水水位低于右側(cè)淡水水位，根據(jù)式(14)得到

(27)

(28)

當(dāng)滿足

(29)

根據(jù)絕對值不等式

(30)

3 數(shù)值算例

3.1 正向滲流算例

假設(shè)存在一個滲流槽，用于模擬海水入侵，槽長L=1.5 m，寬B=0.1 m，高0.6 m，內(nèi)部填充滿含水介質(zhì)，K=0.000 1 m/s，αL=0.2 m，左側(cè)為定水頭和定濃度邊界，左側(cè)海水水位為0.5 m，左側(cè)海水初始 離子濃度為18 000 mg/L；右側(cè)為定水頭邊界，右側(cè)淡水水位為0.35 m，且邊界衰減因子為β，此情況屬于正向滲流類型，其示意圖見圖3。取Δx=0.03 m，Δt=100 s，求經(jīng)過一定迭代次數(shù)后，在滲流槽中海水鹽分?jǐn)U散后穩(wěn)定的濃度與擴(kuò)散距離的關(guān)系。初始及邊界條件為

圖3 正向滲流類型下滲流槽尺寸及邊界條件示意圖

圖4 正向滲流類型下最終濃度隨滲流距離的變化情況

采用有限差分法計算上述算例，在滿足收斂條件的情況下，經(jīng)過6 000步迭代次數(shù)后，得到節(jié)點位置的濃度基本上不隨迭代次數(shù)的增加而改變，則繪制最終濃度隨滲流距離變化的曲線圖，并比較邊界衰減因子β=0與β=1的情況，如圖4所示。 在正向滲流的情況下，右邊界處濃度的衰減對最終穩(wěn)定狀態(tài)的影響是比較大的，不能忽略。

3.2 逆向滲流算例

當(dāng)左側(cè)海水水位為0.35 m、右側(cè)淡水水位為0.5 m時，其他條件均與上述情況相同，此情況為逆向滲流類型，其示意圖見圖5。采用有限差分法計算上述算例，繪制最終濃度隨滲流距離變化曲線圖，并比較邊界衰減因子β=0與β=1的情況，如圖6所示。

圖5 逆向滲流類型下滲流槽尺寸及邊界條件示意圖

圖6 逆向滲流類型下最終濃度隨滲流距離的變化情況

在逆向滲流的情況下，右邊界處濃度的衰減對最終穩(wěn)定狀態(tài)的影響可忽略不計。

4 結(jié)論

利用顯式有限差分法求解一維對流彌散方程的數(shù)值解，并結(jié)合2種不同滲流方向的類型，給出了相應(yīng)的顯式差分式，推導(dǎo)了各自應(yīng)滿足的收斂條件，且通過實際案例進(jìn)行算法演示，求解了濃度擴(kuò)散數(shù)值，最后探討了在不滿足收斂條件的情況下的數(shù)值特性。計算結(jié)果表明:

(1)對于2種不同的類型，其有限差分的迭代關(guān)系式是不同的，要根據(jù)海水入侵的速度方向確定屬于什么類型，選擇相應(yīng)的迭代關(guān)系式，且要預(yù)先計算出流速場，流速場采用本文介紹的計算式算出，然后試代Δx和Δt的值，判定是否滿足收斂條件，根據(jù)收斂條件選擇適當(dāng)?shù)臅r間步長Δt和空間步長Δx。Δx和Δt的值越小則截斷誤差越小，計算的結(jié)果精確度越高，但相應(yīng)的計算量和迭代次數(shù)會增加，在滿足收斂條件下，適當(dāng)?shù)靥岣擀的值會大大減少迭代次數(shù)。采用顯式差分法求解一維對流彌散問題，經(jīng)過一定的迭代次數(shù)后，發(fā)現(xiàn)數(shù)值穩(wěn)定了，或者誤差在允許范圍內(nèi)，所得到的最終結(jié)果便是經(jīng)過相應(yīng)時間擴(kuò)散后的濃度，只要滿足相應(yīng)的收斂條件，其計算結(jié)果是收斂且誤差較低的。

(2)在正向滲流的類型下，右邊界處濃度的衰減對最終穩(wěn)定狀態(tài)的影響是比較大的，不能忽略。在逆向滲流的類型下，右邊界處濃度的衰減對最終穩(wěn)定狀態(tài)的影響可忽略不計。

(3)本顯式有限差分法能夠?qū)崿F(xiàn)對一維海水入侵對流彌散問題的求解，以往對一維對流彌散方程的求解并沒有討論正向和逆向滲流的情況，討論在2種不同滲流方向下發(fā)生鹽分遷移，利用本文所述方法能夠計算出各自穩(wěn)定的濃度情況，原理簡單易操作，且能夠方便了解每一迭代步內(nèi)的濃度空間分布情況，更直觀研究每個迭代步內(nèi)的濃度變化規(guī)律，是求解海水入侵問題的有效方法。