基于Laplace變換的潛水非穩定流模型解的研究

劉 能 勝

(湖北水利水電職業技術學院，武漢 430070)

0 引 言

河渠水位的變化是影響兩岸地下水動態的重要因素，對河渠間潛水的非穩定運動規律的研究，對區域地下水資源評價、灌排渠系設計、地下水動態預報都有著重要的意義[1-4]。其中，地表水體附近潛水非穩定流與滲流規律研究[5-11]已在渠道附近地下水動態過程研究[12-15]中得到廣泛應用。

在之前的研究中，假設河渠完全切割潛水含水層，將河渠視為第一類邊界的河渠—半無限潛水非穩定流模型，大體分為兩類：①河渠水位變化瞬時完成并保持穩定[1-4]，對于這類簡單的邊界條件函數，可利用常用的積分變換方法，求出模型的理論解；②河渠水位持續變化[1,3,5,6]，對于這類問題多將水位的變化過程概化為一常規函數，如線性變化，當函數形式較為簡單時，也可以通過常用的積分變化方法進行模型求解。對于上述兩類問題現已有較為深入的研究，但實際情況中，河渠水位多難以一已知函數形式進行概化，即便某些特定條件下可給出變化過程的具體函數形式，但當函數過于復雜時[10]，其常用的積分變換方法便難以進行模型求解。

后來的研究中為解決復雜函數邊界的河渠—半無限潛水非穩定流模型問題，吳丹[16]等利用傅里葉變換方法對所概化的復雜函數邊界進行處理，并給出了該邊界條件下模型解，但其求解過程過于復雜，且并未驗證解的正確性，而且在實際情況中，其河渠水位變化過程多為一未知函數，故對未知邊界函數形式條件下的河渠—半無限潛水非穩定流模型問題進行相關研究具有實際意義。

針對上述情況，本文采用不考慮邊界條件形式的Laplace變換方法求解，利用Laplace變換中的有關性質，給出模型的解，并與數值解進行對比驗證；同時為了進一步討論該理論解在實際中的運用，本文對河渠水位變化過程采用Lagrange插值進行離散化處理，并借助文獻[16,17]中的算例數據進行驗算，對其應用展開介紹。

1 經典模型及其解

為研究河渠附近潛水非穩定運動，對河渠及附近地段的水文地質條件作如下假設：①潛水含水層均質各向同性，隔水底板水平、平面上無限延伸；②河渠順直且在剖面上完整切割含水層，河渠水位變動過程為時變函數H(t)；③潛水初始水位h(x,0)與河渠水位水平、潛水水流可視為一維流；④忽略上部入滲量，不考慮垂向交換作用。

上述水文地質概念模型，可寫成數學模型(Ⅰ)。

圖1 河渠附近潛水的非穩定運動Fig.1 Transient phreatic motion near rivers and canals

式中：μ為潛水含水層給水度；H為潛水水位，m；K為含水層滲透系數，m/d；H(t)為河渠水位隨時間變化函數，m。

在經典模型中，H(t)=H0，H0為常數，即“河流水位在迅速升高H0后，保持不變”。當含水層厚度較大時，水位變化h(x,t)-h(x,0)≤0.1hm(潛水流的平均厚度，即研究時段始、末潛水流厚度的算術平均值)時，根據Boussinesq方程第一線性化方法，模型(Ⅰ)可改寫成(Ⅱ)。

(Ⅱ)中：u(x,t)=h(x,t)-h(x,0)，a=khm/μ，a，(m2/d)為潛水含水層的導壓系數。

當H(t)=H0，利用Laplace變換，可獲得上述模型的解：

(9)

2 理論模型

2.1 理論模型的求解

實際中的河渠水位變化過程H(t)，雖在小范圍流域內有一定的規律，但由于外界干擾因素過多，其變化過程更具隨機性，再加灌區人工干擾因素，基本無法對水位變化過程進行統一的具體函數表達，在此條件下，本文在不追求H(t)的具體函數形式基礎上，建立不依賴邊界函數的變換過程的Laplace變換方法，求解這類模型的理論解。

對于模型(Ⅲ)，當式(11)過于復雜或未知時，難以采用上述方法直接求解，對此，針對式(11)，不直接求H(t)的象函數，而只將其以算符的形式運用于計算過程中，再利用變換的相關性質和相應的特征函數，進行求解。詳解過程如下。

(Ⅲ)中(10)式的通解為：

(13)

式(13)中，λ1、λ2為待定常數；由邊界條件(11)和(12)，模型(Ⅲ)的特定解為：

(14)

(15)

式中：*為卷積算符。

參照常用Laplace變換中的函數：

對式(16)，根據Laplace變換的“微分定理”，有：

(17)

其中：

結合式(15)、(17)，有：

(18)

對式(18)，將卷積打開，注意h(x,t)和u(x,t)的關系可得：

式(19)即為不考慮H(t)的具體函數形式下模型(Ⅱ)的解；在實際問題的運用中，還需將具體的H(t)代入到式中，進一步換算，才能得到實際問題的解。

對于式(19)，其河渠水位變化過程以H(t)進行概化，而之前的研究中多以一已知常規函數進行處理，但在實際中，河渠水位的變化多難以函數形式進行概化，即便存在，也多為一些非常規函數，給求解帶來了極大的困難，而本文通過對Laplace變換性質的運用，避免了邊界函數的變換過程的復雜性，進一步推動了經典模型的解在實踐中的運用。同時需要指出的是，上述求解過程中，邊界條件H(t)雖未進行形式上的變換，但其實際參與了Laplace變換與逆變換的過程，所以，H(t)應滿足Laplace變換要求，在實際中，河渠水位變化基本是可以滿足變化要求的。

2.2 解析解與數值解的對比分析

為了進一步探討解析解的正確性以及與模型數值解的區別，本文以下述應用中文獻[17]中的實例為背景，利用基于有限差分的GMS軟件[18]對河渠附近地下水進行數值模擬，將河渠作為第一類邊界處理，垂向上為與解析解概化一致，只設置一層，并對潛水一維穩定流的數值解與解析解的計算結果進行比較，計算模型概化中，河渠水位變化過程如表2所述，各參數取值為給水度μ=0.04，滲透系數k=4.26 m/d，計算不同時刻河渠附近地下水位情況。

圖2為河渠水位變化24 h末的數值解和解析解的比較。

圖2 24 h末河渠附近地下水位Fig.2 24 h-end Water Table Near Rivers and Canals

從圖2可以看出，至24 h末，河渠附近地下水水位的數值解和解析解存在一定誤差，經分析，解析解與數值解的偏差率為0.03%，平均相對誤差為0.4%，通過分析解析解求參運用及與數值解的對比可得：

(1)當實際水文地質條件較為簡單時，在合理概化的基礎之上可以直接采用解析法來求取水文地質參數并預測地下水動態。

(2)解析解在野外實際運用時，可以提供一個簡潔的估算方法，也可以利用其所得結果進一步指導勘探和野外試驗工作。

(3)對本文所述的理論模型來說，解析解是精確解，可以用來檢驗數值解法的正確性和比較不同數值解法的精確性和有效性。

同時，解析法的局限性也非常明顯，它不能適用于含水層的各種復雜條件，因此將它利用于實際計算時不可能得到精確的結果，但作為一種進行理論研究和初步估算的工具來說，具有重要的意義。

3 運 用

3.1 Lagrange插值

實際中，復雜的河渠水位變化過程H(t)，顯然難以給出統一的數學表達式；在之前的研究中，多將河渠水位變化過程多概化為階梯狀分段函數，但在實際情況中河道水位的變化應是連續的，階梯狀的河道水位變化過程難以出現，因此，為使河道水位變化過程的概化的實際意義更加合理，本文采用Lagrange插值方法[16,19]，對H(t)進行離散化處理；

對于時段Ti=ti-ti-1，有：

(20)

則有：

(21)

δ(t-ti-1)系Heaviside 函數，具有如下性質：當t

3.2 求 參

在吳丹、陶月贊[16,17]等有關河渠附近潛水非穩定流的研究中，提出了“配線法”和“拐點法”的參數求解方法，推廣了解析解在實際中的運用，本文將基于上述兩種方法，結合相關實例探討本文所求解析解在實際中的運用。

設潛水位變動速度φ(x,t)=?h(x,t)/?t，注意u(x,t)=h(x,t)-h(x,0)，結合卷積定理的微分特性，由(18)式：

(22)

結合式(21)，注意，當t=0時，H(t)=H0及Heaviside 函數的性質，得：

(23)

3.2.1 配線法

對于式(23)，當H0=0∩λ≠0時：

(24)

需要注意的是，配線法在實際運用時，要求河渠初始水位基本保持穩定，這樣從模型概化的角度出發即可認為H0=0。

3.2.2 拐點法

對于式(22)，當H0=0∩λ≠0時，對φ(x,t)關于t再進行一次求導，并令?φ(x,t)/?t=0，求出相應拐點；

(25)

令拐點處的時間為tk，有：

(26)

根據距邊界x處的地下水水位動態，繪出?φ(x,t)/?t～t曲線(此時，x是一確定值)；利用曲線的拐點tk，以及河渠水位實際變動數據而獲得的H0、λ；由(26)式，可計算出含水層參數a，需要注意的是，拐點法應用時要求河渠初始水位要有個迅速的變化過程，這樣從模型概化的角度出發即可認為H0≠0。

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

河渠與潛水之間的補排關系，由河渠水位與潛水水位關系所決定，這也表明，河渠附近潛水非穩定滲流模型，可以用來計算河渠與潛水之間的水量交換；其中，河渠與潛水之間的交換量，可以用單位河渠長度上的交換量表達，即單寬流量q(t)；t時刻的交換強度q(t)，由達西定理：

(27)

如前文，在實際問題的運用中，還需將具體的H(t)帶入到式中，進一步換算，才能得到河渠與潛水之間的水量交換計算公式。若H(t)采用Lagrange插值進行離散化處理，根據式(27)可得：

(28)

式(28)即為河渠與潛水之間在t時刻的交換強度算式；q(t) >0，表示潛水接受河渠補給；q(t)<0，表示潛水排泄河渠。

4 實 例

本文以文獻[16,17]中算例數據進行參數計算。根據文獻所述，安徽淮北平原中部的蒙城縣境內茨淮新河灌區潛水含水層廣泛發育，厚度8 m左右，其下發育有不完全連續黏性土隔水層；潛水位埋深多較淺；且農業區內灌渠系統較完善，灌渠之間的間距多為2 km左右，干渠渠首大部分建設有節制閘；區內有一口國家級地下水位自記觀測井，距離干渠距離約為60 m，觀測井附近地面標高31.02 m。

據文獻[16]表述，2014年7月下旬，由于未處于運行期，干渠水位基本長時間保持不變，符合“配線法”中要求的H0=0；且該時段為無雨期，計算中忽略垂向交換作用；可采用配線法進行求參；數據摘錄如表1。

據文獻[17]表述，2013年8月26日15時至27日15時，干渠節制閘關閘蓄水；其水位變化過程可分為兩個過程，關閘15 min內，水位迅速上升2.00 m，隨后，水位上升速度減緩，至27日15時升幅為0.21 m。觀測孔水位摘錄如表2。該時段，H0=2.0 m、λ=0.21 m/d，滿足拐點法要求的求參條件，可采用拐點法求參數。

由圖3可見，配線法中，實測散點數據基本處與a值為700～900 m2/d，因此a值可取800 m2/d；對于拐點法，其潛水位達到最大變動速度時所對應的時間十分靠近18 h，故拐點取tk=17.2 h，按式(26)求得a=851 m2/d。

表1 潛水水位動態數據(2014.7.25-2014.7.30)與“配線法”計算過程Tab.1 The dynamic data of water table (2014.7.25-2014.7.30) and calculation process of “Curve Fitting Method”

表2 潛水水位動態數據(2013.8.26-2014.8.27)與“拐點法”計算過程Tab.2 The dynamic data of water table (2013.8.26-2014.8.27) and calculation process of “Inflection Point Method”

圖3 配線法及拐點法求a(x=60 m)Fig.3 Calculate a(x=60 m) by means of curve fitting and inflection point methods

根據文獻[16,17]的研究成果，采用配線法及拐點法求得參數a值分別為890 m2/d和860 m2/d，與本文所求結果基本一致，但仍存在一定差異，其主要原因為：

(1)上述兩方法在最后配線及拐點位置確定時均為人為估計，難免會產生一定誤差；

(2)本文進行參數求解與上述文獻進行參數求解時所依據的理論模型及解析解有所區別，具體表現在對河渠水位邊界形式概化上。

5 結 論

本文通過采用不考慮邊界條件形式的Laplace變換方法對未知函數邊界條件下的河渠—半無限潛水非穩定流模型進行求解，并與數值解進行對比驗證，同時展開了相關應用的探討，得出以下結論：

(1)針對半無限域河渠附近潛水非穩定運動經典模型中河渠水位邊界條件概化的局限性，本文采用不考慮邊界條件形式的Laplace變換方法求解，可得到形式較為簡單的解析解，此外，Laplace變換在模型求解中可避免邊界條件函數變化所帶來的象函數求逆的困難，對于未知邊界函數的模型求解有著重要的意義。

(2)河渠水位變化過程往往比較復雜時，利用Lagrange插值的方法對河渠水位過程進行離散化處理，同時利用配線法及拐點法進行參數計算，參考相關文獻，結果基本一致，方法可行，此外，插值函數的基礎上，推導了河渠與潛水之間交換量計算，可為實際問題求解提供更簡便的途徑。

(3)對比數值解和解析解可以發現，當實際水文地質條件較為簡單時，在合理概化的基礎之上可以直接采用解析法來求取水文地質參數并預測地下水動態，但需要注意的是解析法的局限性也非常明顯，它不能適用于含水層的各種復雜條件，此外模型是建立在一維概化的條件下，因此將它利用于實際計算時不可能得到精確的結果，但作為一種進行理論研究和初步估算的工具來說，具有重要的意義。

(4)為了計算的便捷性，本文對水位變化過程只進行了簡單的插值離散處理，隨著計算機技術的發展，可以對長系列的水文數據進行函數擬合，考慮到水位變化在實際情況中還存在其他情形，需進一步研究各情況下理論解的實際應用。

(5)在河渠—半無限潛水非穩定流模型中，除河渠水位變換對附近潛水具有重要影響外，其垂向水量交換及含水層非均質特性也是重要的影響因素，筆者認為，在模型中可將垂向水量交換概化為源匯項，假設其強度為一穩定值；此外對于含水層非均質特性問題，可以先著眼于層狀非均質模型的研究，受限于文章篇幅，這里不再作詳細展開，但應作為后續研究的重點。