基于Laplace變換的潛水非穩(wěn)定流模型解的研究

劉 能 勝

(湖北水利水電職業(yè)技術(shù)學(xué)院，武漢 430070)

0 引 言

河渠水位的變化是影響兩岸地下水動(dòng)態(tài)的重要因素，對(duì)河渠間潛水的非穩(wěn)定運(yùn)動(dòng)規(guī)律的研究，對(duì)區(qū)域地下水資源評(píng)價(jià)、灌排渠系設(shè)計(jì)、地下水動(dòng)態(tài)預(yù)報(bào)都有著重要的意義[1-4]。其中，地表水體附近潛水非穩(wěn)定流與滲流規(guī)律研究[5-11]已在渠道附近地下水動(dòng)態(tài)過(guò)程研究[12-15]中得到廣泛應(yīng)用。

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

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

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

1 經(jīng)典模型及其解

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

上述水文地質(zhì)概念模型，可寫成數(shù)學(xué)模型(Ⅰ)。

圖1 河渠附近潛水的非穩(wěn)定運(yùn)動(dòng)Fig.1 Transient phreatic motion near rivers and canals

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

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

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

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

(9)

2 理論模型

2.1 理論模型的求解

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

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

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

(13)

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

(14)

(15)

式中：*為卷積算符。

參照常用Laplace變換中的函數(shù)：

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

(17)

其中：

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

(18)

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

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

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

2.2 解析解與數(shù)值解的對(duì)比分析

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

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

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

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

(1)當(dāng)實(shí)際水文地質(zhì)條件較為簡(jiǎn)單時(shí)，在合理概化的基礎(chǔ)之上可以直接采用解析法來(lái)求取水文地質(zhì)參數(shù)并預(yù)測(cè)地下水動(dòng)態(tài)。

(2)解析解在野外實(shí)際運(yùn)用時(shí)，可以提供一個(gè)簡(jiǎn)潔的估算方法，也可以利用其所得結(jié)果進(jìn)一步指導(dǎo)勘探和野外試驗(yàn)工作。

(3)對(duì)本文所述的理論模型來(lái)說(shuō)，解析解是精確解，可以用來(lái)檢驗(yàn)數(shù)值解法的正確性和比較不同數(shù)值解法的精確性和有效性。

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

3 運(yùn) 用

3.1 Lagrange插值

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

對(duì)于時(shí)段Ti=ti-ti-1，有：

(20)

則有：

(21)

δ(t-ti-1)系Heaviside 函數(shù)，具有如下性質(zhì)：當(dāng)t

3.2 求 參

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

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

(22)

結(jié)合式(21)，注意，當(dāng)t=0時(shí)，H(t)=H0及Heaviside 函數(shù)的性質(zhì)，得：

(23)

3.2.1 配線法

對(duì)于式(23)，當(dāng)H0=0∩λ≠0時(shí)：

(24)

需要注意的是，配線法在實(shí)際運(yùn)用時(shí)，要求河渠初始水位基本保持穩(wěn)定，這樣從模型概化的角度出發(fā)即可認(rèn)為H0=0。

3.2.2 拐點(diǎn)法

對(duì)于式(22)，當(dāng)H0=0∩λ≠0時(shí)，對(duì)φ(x,t)關(guān)于t再進(jìn)行一次求導(dǎo)，并令?φ(x,t)/?t=0，求出相應(yīng)拐點(diǎn)；

(25)

令拐點(diǎn)處的時(shí)間為tk，有：

(26)

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

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

河渠與潛水之間的補(bǔ)排關(guān)系，由河渠水位與潛水水位關(guān)系所決定，這也表明，河渠附近潛水非穩(wěn)定滲流模型，可以用來(lái)計(jì)算河渠與潛水之間的水量交換；其中，河渠與潛水之間的交換量，可以用單位河渠長(zhǎng)度上的交換量表達(dá)，即單寬流量q(t)；t時(shí)刻的交換強(qiáng)度q(t)，由達(dá)西定理：

(27)

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

(28)

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

4 實(shí) 例

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

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

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

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

表1 潛水水位動(dòng)態(tài)數(shù)據(jù)(2014.7.25-2014.7.30)與“配線法”計(jì)算過(guò)程Tab.1 The dynamic data of water table (2014.7.25-2014.7.30) and calculation process of “Curve Fitting Method”

表2 潛水水位動(dòng)態(tài)數(shù)據(jù)(2013.8.26-2014.8.27)與“拐點(diǎn)法”計(jì)算過(guò)程Tab.2 The dynamic data of water table (2013.8.26-2014.8.27) and calculation process of “Inflection Point Method”

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

根據(jù)文獻(xiàn)[16,17]的研究成果，采用配線法及拐點(diǎn)法求得參數(shù)a值分別為890 m2/d和860 m2/d，與本文所求結(jié)果基本一致，但仍存在一定差異，其主要原因?yàn)椋?/p>

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

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

5 結(jié) 論

本文通過(guò)采用不考慮邊界條件形式的Laplace變換方法對(duì)未知函數(shù)邊界條件下的河渠—半無(wú)限潛水非穩(wěn)定流模型進(jìn)行求解，并與數(shù)值解進(jìn)行對(duì)比驗(yàn)證，同時(shí)展開了相關(guān)應(yīng)用的探討，得出以下結(jié)論：

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

(2)河渠水位變化過(guò)程往往比較復(fù)雜時(shí)，利用Lagrange插值的方法對(duì)河渠水位過(guò)程進(jìn)行離散化處理，同時(shí)利用配線法及拐點(diǎn)法進(jìn)行參數(shù)計(jì)算，參考相關(guān)文獻(xiàn)，結(jié)果基本一致，方法可行，此外，插值函數(shù)的基礎(chǔ)上，推導(dǎo)了河渠與潛水之間交換量計(jì)算，可為實(shí)際問(wèn)題求解提供更簡(jiǎn)便的途徑。

(3)對(duì)比數(shù)值解和解析解可以發(fā)現(xiàn)，當(dāng)實(shí)際水文地質(zhì)條件較為簡(jiǎn)單時(shí)，在合理概化的基礎(chǔ)之上可以直接采用解析法來(lái)求取水文地質(zhì)參數(shù)并預(yù)測(cè)地下水動(dòng)態(tài)，但需要注意的是解析法的局限性也非常明顯，它不能適用于含水層的各種復(fù)雜條件，此外模型是建立在一維概化的條件下，因此將它利用于實(shí)際計(jì)算時(shí)不可能得到精確的結(jié)果，但作為一種進(jìn)行理論研究和初步估算的工具來(lái)說(shuō)，具有重要的意義。

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

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