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

1 研究背景

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

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

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

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

2 研究方法

2.1 一維對(duì)流彌散方程

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

(1)

(2)

Dl=αL|Vx|+D*

(3)

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

Dl=αL|Vx|

(4)

2.2 劃分網(wǎng)格

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

2.3 微分方程的顯式差分化

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

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

(5)

(6)

(7)

(8)

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

(9)

(10)

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

(11)

令

(12)

則可化簡(jiǎn)為

(13)

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

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

(14)

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

(15)

(16)

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

(17)

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

2.4 差分格式的收斂性判定

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

記

(18)

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

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

(19)

令

(20)

有

(21)

當(dāng)滿足

(22)

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

(23)

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

(24)

如此類推，對(duì)于某個(gè)時(shí)刻tM有

(25)

(26)

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

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

(27)

(28)

當(dāng)滿足

(29)

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

(30)

3 數(shù)值算例

3.1 正向滲流算例

假設(shè)存在一個(gè)滲流槽，用于模擬海水入侵，槽長(zhǎng)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ì)算上述算例，在滿足收斂條件的情況下，經(jīng)過6 000步迭代次數(shù)后，得到節(jié)點(diǎn)位置的濃度基本上不隨迭代次數(shù)的增加而改變，則繪制最終濃度隨滲流距離變化的曲線圖，并比較邊界衰減因子β=0與β=1的情況，如圖4所示。 在正向滲流的情況下，右邊界處濃度的衰減對(duì)最終穩(wěn)定狀態(tài)的影響是比較大的，不能忽略。

3.2 逆向滲流算例

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

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

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

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

4 結(jié)論

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

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

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

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