復(fù)雜函數(shù)邊界控制下的潛水非穩(wěn)定流模型及解的應(yīng)用

吳　丹，陶月贊，林　飛

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

1　研究背景

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

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

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

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

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

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

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

圖1　河渠附近潛水滲流場(chǎng)

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

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

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

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

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

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

當(dāng)h（x，t）-h（x，0）≤0.1hm（hm為計(jì)算期間的潛水流平均厚度，實(shí)際中的潛水問題多能符合該條件）時(shí)，可利用 Boussinesq方程的第一線性化方法；同時(shí)，令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ù)時(shí)，即為J.G.Ferris模型；此時(shí)模型的解為［1-3］：

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

3　模型理論解

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

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

所以，由式（5），有：

上式通解為：

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

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

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

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

式中*為卷積算符。

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

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

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

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

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

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

根據(jù)Lagrange線性插值原理，在時(shí)段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ù)對(duì)應(yīng)的解將式（17）帶入式（15），有：

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

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

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

5　解的應(yīng)用

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

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

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

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

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

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

圖2　適線法求a

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

5.2河渠與潛水之間交換量計(jì)算

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

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

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

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

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

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

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

6　實(shí)例研究

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

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

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

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

7　結(jié)語

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

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

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

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

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

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