吳　丹，陶月贊，林　飛

（合肥工業(yè)大學(xué) 土木與水利工程學(xué)院，安徽 合肥　230009）

1　研究背景

將完全切割潛水含水層的河渠水位作為第一類邊界，邊界一側(cè)潛水非穩(wěn)定滲流問題，是地下水滲流力學(xué)中的最經(jīng)典問題之一［1-3］。該問題也屬帶自由面滲流問題，所以一直是地下水滲流力學(xué)領(lǐng)域中的研究難點與熱點之一［4］。上述問題的解，可為灌溉渠系設(shè)計提供重要的計算依據(jù)［1］：通過調(diào)控渠系水位將影響區(qū)地下水水位控制在合理區(qū)間，是現(xiàn)代許多灌區(qū)的基本要求［5-6］。

就邊界條件而言，河渠水位變化過程不僅受河渠自身水力特征和區(qū)域水文氣象條件控制，而且受到明顯的人工調(diào)控影響，具強烈的隨機性［7-8］；對這種復(fù)雜多變的時變函數(shù)所構(gòu)成的邊界條件，隨著對其概化與處理的方法不同，模型的求解方法以及所獲得的解也有所差別。

這類問題的經(jīng)典模型——J.G.Ferris模型（1950年），假定河渠水位變化瞬時完成之后長期保持穩(wěn)定［1，3，9-10］，即下文模型（Ⅰ）中的邊界條件H（t）為常數(shù)H0；對河渠水位復(fù)雜多變條件下的模型，利用常用積分變換方法求解時，按對邊界處理方式與求解方法不同，現(xiàn)行方法主要可分為三大類：

方法一。H（t）是如時間線性函數(shù)等相對簡單的函數(shù)，可通過常用的積分變換方法對模型直接求解［1，3，9-10］；該方法對一些特定情形（如經(jīng)典模型所假設(shè)的情形）有較好的適用性，但不具普遍應(yīng)用意義。

方法二。H（t）采用水位或水位變速“分段等值”方法離散化處理，依據(jù)經(jīng)典模型的解，利用疊加原理構(gòu)建問題的解［1］。這類方法就求解過程而言，方法比較簡便；但當(dāng)H（t）比較復(fù)雜時，劃分的時段數(shù)較多，計算偏繁瑣。

方法三。H（t）形式復(fù)雜，包括一些初始條件較復(fù)雜問題，在求解過程中或解的表達式中，需借助不常用的特殊函數(shù)或是構(gòu)建專用的特殊函數(shù)［4-6，11-14］，這為應(yīng)用帶來不便。綜合上述方法優(yōu)點，針對這類問題的數(shù)學(xué)模型特點，在不依賴邊界函數(shù)的變換過程的條件下，利用Fourier變換并依據(jù)卷積定義和卷積的微分性質(zhì)，給出模型的理論解；對實際河渠水位過程，進行Lagrange線性插值離散，將插值函數(shù)代入理論解，給出關(guān)于問題的實際解；這一求解過程比較簡單，所獲得的問題實際解，不僅形式相對簡單，而且是由常用函數(shù)表達。

2　經(jīng)典模型及基本解

圖1　河渠附近潛水滲流場

河渠及附近地段的水文地質(zhì)條件（圖1），可概括如下：

①具水平隔水底板的潛水含水層均質(zhì)各向同性、平面無限延展；

②河渠水位變動過程，為時變函數(shù)H（t）；

③初始潛水位h（x，0）與河渠水位水平；

④河渠順直且在剖面上基本完整切割含水層，潛水水流可視為一維流。該水文地質(zhì)概念模型所對應(yīng)的數(shù)學(xué)模型，可寫成模型（Ⅰ）：

式中：μ為含水層給水度；h為潛水水位，m；K為滲透系數(shù)，m/d；H（t）為第一類邊界（河渠水位）隨時間變化的函數(shù)，m。

當(dāng)h（x，t）-h（x，0）≤0.1hm（hm為計算期間的潛水流平均厚度，實際中的潛水問題多能符合該條件）時，可利用 Boussinesq方程的第一線性化方法；同時，令u（x，t）=h（x，t）-h（x，0），a=khm/μ，a為潛水含水層的導(dǎo)壓系數(shù)，m2/d；可將模型（Ⅰ）轉(zhuǎn)化為模型（Ⅱ）。

當(dāng)H（t）=H0，H0為已知常數(shù)時，即為J.G.Ferris模型；此時模型的解為［1-3］：

式中，erfc（z）是余誤差函數(shù)，z=x/2（at）1/2。式（9）也是這類問題的基本解。

3　模型理論解

對于模型（Ⅱ），在不依賴H（t）具體變換過程的條件下，建立這類模型的理論解。

由邊界條件與變量的變化范圍，根據(jù)Fourier變換的性質(zhì)與特點，對模型（Ⅰ）求關(guān)于x的Fourier正弦變換。為u關(guān)于x的Fourier變換，ω、F和F-1分別為變換算子、變換與逆變換算符。有：

所以，由式（5），有：

上式通解為：

對式（10），求正弦逆變換，注意正弦與余弦逆變換關(guān)系，并適時交換積分次序：

式（12）是在H（t）還未明確條件下模型（Ⅱ）的解；也即對任意H（t），式（12）是模型的理論解。具體到實際問題應(yīng)用，還需要將已知的H（t）代入，再進一步展開，才能獲得實際問題的解。

為獲得更為便捷的求解方法，并結(jié)合本領(lǐng)域的應(yīng)用習(xí)慣，以下依據(jù)積分變換性質(zhì)及卷積定理，給出以概率密度函數(shù)表達的解。

據(jù)卷積定義，由式（12），有：

式中*為卷積算符。

由卷積的微分性質(zhì)：

式（15）實質(zhì)上也是模型（Ⅱ）的另外一種形式的理論解。

上述求解過程中，邊界條件H（t）一直未進行形式上的變換，也即在H（t）屬未知的情況下進行的；這一解法可為H（t）復(fù)雜條件下問題求解提供支持。

需要說明的是，上述方法求解過程中，雖然H（t）未進行形式上的變換，但H（t）實質(zhì)上參與了變換與逆變換過程；所以H（t）必須滿足Fourier變換要求，即：函數(shù)在任意區(qū)間滿足Dirichlet條件，在無限區(qū)間上絕對可積；就河渠水位實際變化過程而言，通常都可滿足上述要求的。

4　河渠水位過程離散及其解

4.1Lagrange線性插值實際中，復(fù)雜的水位變化過程H（t）顯然難以給出統(tǒng)一的數(shù)學(xué)表達式；此條件下，不依賴H（t）的具體函數(shù)形式，只根據(jù)H（t）的實測過程，采用Lagrange線性插值方法對其進行離散化處理。

根據(jù)Lagrange線性插值原理，在時段Ti=ti-ti-1，有：

則有：

ε（t-ti-1）系Heaviside 函數(shù)［12］， 具有如下性質(zhì)：當(dāng)t＜ti-1、ε（t-ti-1）=0，當(dāng)t≥ti-1、ε（t-ti-1）=1。

4.2線性插值函數(shù)對應(yīng)的解將式（17）帶入式（15），有：

注意ε（t）函數(shù)的性質(zhì)及u（x，t）=h（x，t）-h（x，0），整理式（18），可得：

式（19）是在河渠邊界控制下的半無限潛水含水層中，對水位變化過程H（t）采用Lagrange線性插值離散化的條件下，潛水非穩(wěn)定滲流過程的解析解。

當(dāng)然，也可將式（17）帶入式（13），再通過分部積分等步驟進一步展開，以獲得最終形式的解；但由式（15）而獲得式（19），求解過程顯然要簡明得多。

5　解的應(yīng)用

河渠潛水滲流模型研究的重要目的之一是以其為工具，利用潛水位變動數(shù)據(jù)計算模型參數(shù)a。

5.1模型參數(shù)計算潛水位變動速度φ（x，t）=?h()x，t ?t，由式（16）：

以下僅討論H0、λ1的情形。

（1）當(dāng)H0=0∩λ1≠0的配線法。由式（20）：

對距離邊界為x的觀測孔（x為確定值），z=x/2（at）1/2，首先建立不同a值的erf（z）～t理論曲線圖族；由實測地下水水位，制出φ（x，t）～t的曲線；由式（19），當(dāng)φ（x，t）～t所依賴的含水層參數(shù)a值，與erf（z）～t理論曲線圖族中某條曲線a值相等時，則，兩條曲線形態(tài)完全相同，僅相差一個常數(shù)λ1倍，也即兩條曲線應(yīng)該完全重合。

因此，根據(jù)地下水水位實測數(shù)據(jù)，建φ（x，t）～t曲線，將之與erf（z）～t理論曲線族進行適線，就可確定含水層的a值（圖2）。

圖2　適線法求a

（2）當(dāng)H0≠0∩λ1=0。式（20）轉(zhuǎn)為J.G.Ferris模型及其解，解的利用有許多文獻，在此不作贅述。

5.2河渠與潛水之間交換量計算

5.2.1算式河渠與潛水之間的補排關(guān)系，由河渠水位與潛水水位關(guān)系所決定，這也表明，河渠附近潛水非穩(wěn)定滲流模型，可以用來計算河渠與潛水之間的水量交換；僅討論H0、λ1的情形。

河渠與潛水之間的交換量，用單位河渠長度上的交換量表達，即單寬流量Q（t）（m3·d-1·km-1）；t時刻的交換強度q（t），由Darcy定理：

式中：hm是潛水流厚度，a=khm/μ。

式（22）是河渠與潛水之間在t時刻的交換強度算式；q＞0，表示潛水接受河渠補給；q＜0，表示潛水排泄河渠。

5.2.2λ段的影響規(guī)律由式（19），去掉等式右端第2項，轉(zhuǎn)化為J.G.Ferris模型的交換強度算式，對應(yīng)的是河渠水位瞬間變動H0所獲得的水量交換強度，記為qH；式（19）去掉等式右端第1項，對應(yīng)的是H0=0條件下、僅由河渠水位變化段所獲得的交換強度，記為qλ。

λ段的影響，可用qλ/qH來反映：

假設(shè)水位在變動段獲得的水位升幅也為H0（即λt等于H0），由式（23），對交換量的貢獻，λt是H0的2倍；也即水位變動段的水位升幅雖然也為H0，但其對交換量的貢獻，是不考慮變化過程影響而計算出效果的2倍；這是水位在t時段內(nèi)緩慢上升至λt的過程中，被不斷抬升的水位持續(xù)增強補給作用的過程累積效應(yīng)所至。

6　實例研究

安徽淮北平原中部的蒙城縣境內(nèi)茨淮新河灌區(qū)，以粉細砂為主的潛水含水層分布廣泛，厚度8 m左右，底部一般發(fā)育有不完全連續(xù)黏性土層；由于潛水位埋深淺，為2.5～3.0m；區(qū)內(nèi)農(nóng)田灌溉渠系比較完善，干渠基本深切至隔水底板、干渠的渠間距為2 km，干渠渠首多有節(jié)制閘控制；一眼國家級地下水位自記觀測井，距離干渠直線距離為60 m處，觀測井附近地面標高31.02 m。

2014年7月中下旬，干渠長時間處于非運行狀態(tài)，水位基本保持不變；25日至30日，一直處于無雨期，也無灌溉回滲影響，計算中不考慮潛水垂向交換作用；例中，時間段較長，數(shù)據(jù)量也較大，觀測孔水位按12h摘錄，如表1。

表1　潛水水位動態(tài)數(shù)據(jù)(2014.7.25—2014.7.30)與“配線法”計算過程

該時段，H0=0，不符合用拐點法求參數(shù)的條件，但時間段較長，數(shù)據(jù)量足夠大，可依據(jù)這段時間潛水變動速度隨時間變化曲線，采用配線法求算；由上表數(shù)據(jù)，計算出的φ（x，t）～t與erf（z）～t理論曲線圖中的曲線族進行適線（如圖2），適線結(jié)果a值為890 m2/d；這與該地區(qū)相關(guān)文獻［10］的計算結(jié)果（854 m2/d）基本一致。

7　結(jié)語

對水位變化過程復(fù)雜的河渠邊界附近潛水非穩(wěn)定滲流問題，在對水位過程離散化的基礎(chǔ)上，利用Fourier變換給出問題的解，并對解的應(yīng)用進行了相應(yīng)的探討，所得結(jié)論如下：

（1）河渠水位變化過程往往比較復(fù)雜且難以給出具體函數(shù)表達式，利用Fourier變換，在不依賴邊界函數(shù)變換的條件下，獲得理論模型通用解；再采用Lagrange線性插值方法，對實際水位過程離散化處理，將插值函數(shù)代入理論模型通用解，可獲得由比較常用函數(shù)組成的解，且形式較簡單。

（2）在利用Fourier變換求解過程中，可充分利用如卷積定理和卷積的微分性質(zhì)等變換性質(zhì)，由此獲得如文中式（11）的理論通用解，可為實際問題求解提供更簡便的途徑。

（3）河渠水位在其變動過程中，被不斷抬升的河渠水位（以上升為例），使河渠對地下水的補給作用呈持續(xù)增大并伴生有過程累積效應(yīng)；在計算河渠與潛水之間水量交換時，采用水位穩(wěn)定或分段穩(wěn)定的計算方法，將因忽略累積效應(yīng)而導(dǎo)致計算偏差，偏差幅度為2λ1t/H0～λ1t/（H0+λ1t）。

（4）依據(jù)模型解，利用地下水水位動態(tài)監(jiān)測數(shù)據(jù)，可計算含水層參數(shù)a；在H0對地下水水位作用不明顯時段，利用地下水位變動速度隨時間變化過程與理論曲線進行配線，實例表明，方法可行。

值得指出的是，文中的模型及其解，僅適用模型（Ⅰ）所對應(yīng)的水文地質(zhì)條件，即一條直線河渠控制的半無限域潛水滲流問題；隨著水文地質(zhì)條件變化而應(yīng)采用相應(yīng)的滲流問題模型，如有多條平行河渠（如灌區(qū)內(nèi)的支渠）相互干擾時，應(yīng)采用河間地塊潛水滲流問題模型。