楊鷺茜， 蘇 燕， 顧承宇

(1. 福州大學(xué)土木工程學(xué)院， 福建 福州 350108； 2. 臺(tái)灣海洋大學(xué)工學(xué)院， 臺(tái)灣 基隆 20224)

0 引言

Trefftz法是快速求解邊界值問題最常用的邊界無網(wǎng)格法之一， 其應(yīng)用越來越廣泛， 近似解可表示為滿足控制方程函數(shù)的線性組合[1]. 根據(jù)Kita和Kamiya的研究， Trefftz法分為直接法和間接法[2]. Kita等[3]描述Trefftz法在求解三維Poisson方程中的應(yīng)用， 使用笛卡爾坐標(biāo)系中的多項(xiàng)式函數(shù)來近似含有未知函數(shù)的不均勻項(xiàng)， 從而確定泊松方程的解. Trefftz法在工程中的應(yīng)用主要有Laplace方程、 雙諧波方程和二維邊界檢測問題[4-5]. Li等[6-7]對Trefftz法、 配點(diǎn)法和其他邊界方法進(jìn)行全面比較， 總結(jié)出Trefftz法是求解方程最簡便的計(jì)算方法， 該法提供最精確的方程解和最佳的數(shù)值穩(wěn)定性. 為了獲得線性方程組的精確解， Liu[8-11]對 Trefftz法做出改進(jìn)， 將單個(gè)特征長度合并到T型完整基底中， 使得到的線性方程組的條件數(shù)明顯減少. 除單尺度Trefftz法外， Ku等[12]基于標(biāo)量同倫方法提出一般動(dòng)力學(xué)法， 并證明求解線性代數(shù)方程組的數(shù)值穩(wěn)定性， 尤其是對于病態(tài)問題系統(tǒng)， 可得到病態(tài)條件極大的線性方程組在三維 Laplace問題上的解. 當(dāng)前針對沿海地區(qū)地下水的研究較少， 其演化規(guī)律對許多濱海工程具有重要影響. 沿海越流含水層系統(tǒng)的承壓含水層為自由含水層所覆蓋， 且在它們之間存在著垂直方向的弱透水層. Jiao等[13]提出一個(gè)解析解， 用于描述沿海越流含水層地下水受潮汐影響. Li等[14]依據(jù)線性Boussinesq方程提出另一個(gè)近似解析解， 用于研究沿海越流含水層地下水動(dòng)態(tài). Li等[15]將推導(dǎo)的解析解用于研究越流含水層的滲漏和儲(chǔ)水系數(shù)對地下水位變化的影響， 并發(fā)現(xiàn)來自含水層的滲漏會(huì)對承壓含水層的地下水位產(chǎn)生較大的影響.

由于Trefftz法最初是用于處理歐氏空間的邊界值問題， 因此難以找到Trefftz法應(yīng)用于求解與時(shí)間相關(guān)的問題. 本文用Trefftz法首次建立沿海越流含水層系統(tǒng)中承壓地下水流的暫態(tài)模型, 在時(shí)空坐標(biāo)系統(tǒng)的基礎(chǔ)上, 假定時(shí)間為自變量， 時(shí)空坐標(biāo)系中的初始條件和邊界條件均可視為時(shí)空邊界上的邊界條件， 通過對滿足部分邊界配點(diǎn)控制方程的Trefftz基底進(jìn)行迭加， 來逐步近似在時(shí)空坐標(biāo)系中得到的數(shù)值解. 文中還運(yùn)用該方法對一些問題進(jìn)行驗(yàn)證， 包括數(shù)值解與解析解的比較和精度研究， 以及沿海地下水的數(shù)值模擬.

1 Trefftz法基礎(chǔ)理論

1.1 Trefftz基底推導(dǎo)

在均質(zhì)、 各向同性、 等厚的越流含水層系統(tǒng)中， 承壓地下水控制方程如下：

(1)

式中：h為承壓含水層總水頭；x為距離；t為時(shí)間；T為導(dǎo)水系數(shù)；L為越流系數(shù)；S為儲(chǔ)水系數(shù).

為了求得上述控制方程的暫態(tài)解， 應(yīng)用分離變量法做出以下假定

h(x,t)=X(x)K(t)

(2)

將式(2)代入式(1)， 得

(3)

為了簡單易行地推導(dǎo)， 首先引入以下形式

(4)

(5)

將式(4)與式(5)代入式(3)， 得

(6)

為求得該方程的特征值， 采用常數(shù)p保證常數(shù)的正負(fù). 分λ>0，λ=0，λ<0等3種情況.

情況1λ>0， 即λ=p2(0

(7)

因此

X(x)=c1exp(px)+c2exp(-px)

(8)

(9)

將式(8)與式(9)代入式(2)， 得

(10)

K=d2

(11)

將式(9)與式(13)代入式(2)， 得

(12)

式中：c1,c2,d1,d2,A1,A2,B1,B2為任意常數(shù).

情況2λ=0

(13)

因此

X(x)=c3x+c4

(14)

(15)

將式(14)與式(15)代入式(2)， 得

(16)

式中：c3,c4,d3,A3,B3為任意常數(shù).

情況3λ<0， 即λ=-p2(0

(17)

所以

X(x)=c5cos(px)+c6sin(px)

(18)

(19)

將式(18)與式(19)代入式(2)， 得

(20)

式中：c5,c6,d4,d2,A4,B4為任意常數(shù).

綜上三種情況的分析， 利用加法定理對上述所有解進(jìn)行迭加， 得到一個(gè)線性無關(guān)的沿海地下水控制方程的近似解

(21)

式中：A01,B01,A02,B02,A1p,B1p,A2p,B2p為任意常數(shù)；w為適用于控制方程近似解基底函數(shù)的階數(shù).

1.2 時(shí)空配點(diǎn)原理

基于時(shí)空坐標(biāo)系， 本文將時(shí)間定義為自變量， 采用時(shí)空坐標(biāo)系對沿海地區(qū)地下水進(jìn)行暫態(tài)建模. 本文建立的二維坐標(biāo)系， 包括時(shí)間上的一維和空間上的一維， 以求解一維的沿海地區(qū)地下水問題， 時(shí)間和空間邊界處的水頭值均在研究域內(nèi)的時(shí)空邊界上給出, 從而使原始的一維問題轉(zhuǎn)化為二維問題.

從越流含水層系統(tǒng)地下水控制方程式(1)開始， 邊界條件與初始條件如下

圖1 Trefftz法時(shí)空配點(diǎn)示意圖Fig.1 Schematic diagram of Trefftz method with the space-time collocation method

式中：f(x,t)為給定的定水頭邊界條件函數(shù)；v(x, 0)為給定的初始條件函數(shù)； ΓD為定水頭邊界條件所在邊界.

在時(shí)空坐標(biāo)系統(tǒng)中， 將空間軸(x軸)的初始邊界離散為n1點(diǎn)， 將時(shí)間軸(t軸)的邊界離散為n2點(diǎn), 如圖1所示.

根據(jù)不同的邊界條件， 分別將離散后的空間邊界點(diǎn)代入給定的定水頭邊界條件函數(shù)f(x,t)， 空間邊界點(diǎn)代入給定的初始條件函數(shù)v(x, 0)， 將計(jì)算出的邊界值和初始值代入控制方程的近似解表達(dá)式(21)， 并聯(lián)立線性代數(shù)方程組, 以Aα=B的矩陣運(yùn)算形式表示.

(23)

式中： 矩陣A是由Trefftz基底所組成的尺度為aa×bb矩陣， 矩陣B是由邊界值和初始值所組成的尺度為aa×1的矩陣， 待定系數(shù)α為尺度為bb×1的矩陣. 其中aa=2n1+n2，bb=4(w+1).

因?yàn)門refftz基底本身滿足沿海地區(qū)承壓地下水控制方程， 所以只要時(shí)間軸的邊界點(diǎn)和空間軸的初始邊界點(diǎn)滿足各自的邊界條件， 通過Matlab左除運(yùn)算， 由α=Α/Β即可求解待定系數(shù).

將時(shí)空坐標(biāo)系統(tǒng)內(nèi)任意點(diǎn)代入Trefftz基底所構(gòu)成的近似解中， 由Β=Αα即可計(jì)算出任意位置的水頭值h(xi,tj).

2 數(shù)值模擬分析

2.1 解析解數(shù)值驗(yàn)證

本案例通過一維均質(zhì)越流含水層地下水的模型， 以定水頭作為邊界條件進(jìn)行數(shù)值驗(yàn)證， 說明了Trefftz法與時(shí)空坐標(biāo)系統(tǒng)相結(jié)合求解一類邊界條件下承壓地下水問題的可行性.

研究區(qū)域表達(dá)式如下：

Ω={(x,t)|0≤x≤Ln}

(24)

下式是滿足越流含水層地下水控制方程式(1)的解析解， 將其作為研究區(qū)域Ω的第一類邊界條件與初始條件， 以驗(yàn)證數(shù)值模式.

(25)

現(xiàn)設(shè)5 km的承壓含水層的導(dǎo)水系數(shù)T=1.25 m2·h-1， 儲(chǔ)水系數(shù)S=2×10-4， 越流系數(shù)L=0.005 d-1. 利用Trefftz法計(jì)算0～5 h內(nèi)承壓地下水位的變化過程， 并與解析解(25)進(jìn)行比較， 得到二者的最大絕對誤差. 其中， 數(shù)值驗(yàn)證的空間邊界點(diǎn)n1與時(shí)間邊界點(diǎn)n2均取51點(diǎn)，階數(shù)為10階.

圖2顯示了0～5 h范圍內(nèi)地下水問題數(shù)值解和解析解的變化情況， 從圖2可看出， 數(shù)值解與解析解的變化趨勢非常吻合， 有力地說明Trefftz法應(yīng)用于研究越流含水層地下水的可行性. 在不同時(shí)間內(nèi)， 承壓地下水位與距離均呈現(xiàn)線性增長關(guān)系， 且隨著時(shí)間的延長， 地下水位逐步降低至零， 達(dá)到穩(wěn)定狀態(tài)， 這符合實(shí)際地下水流動(dòng)的演變規(guī)律. 為驗(yàn)證Trefftz法求解結(jié)果的準(zhǔn)確性， 將數(shù)值解與解析解進(jìn)行比較， 見圖3， 結(jié)果表明， 二者的最大絕對誤差8.90×10-16， 僅在初始時(shí)刻存在輕微的誤差波動(dòng)， 且隨著時(shí)間和距離的增加， 其誤差始終維持在10-16的高精度范圍內(nèi)， 未出現(xiàn)不收斂的現(xiàn)象.

圖2 數(shù)值解與解析解的變化情況Fig.2 Change in numerical solution and analytical solution

圖3 數(shù)值解與解析解的誤差分析Fig.3 Error analysis of numerical solution and analytic solution

2.2 沿海地下水模擬

本案例假設(shè)某沿海地區(qū)越流含水層系統(tǒng)均質(zhì)、 各向同性， 其自由含水層與承壓含水層的地下水通過中部弱透水層在垂直方向發(fā)生滲漏[13]. 考慮到自由含水層的厚度遠(yuǎn)大于海洋潮汐波動(dòng)的幅度， 合理地忽略自由含水層所受的潮汐波動(dòng)的影響， 并假定其地下水位與平均海平面保持一致， 將平均海平面作為含水層系統(tǒng)的基準(zhǔn)面， 選擇x軸原點(diǎn)位于平均海平面和海灘的交界處， 向內(nèi)陸延伸方向?yàn)檎较颍?見圖4.

圖4 沿海越流含水層示意圖Fig.4 Schematic diagram of coastal leakage aquifers

采用Trefftz法對該沿海越流含水層系統(tǒng)的承壓地下水進(jìn)行波動(dòng)特性模擬測算， 并與有限差分的計(jì)算結(jié)果進(jìn)行趨勢擬合， 驗(yàn)證Trefftz法應(yīng)用于沿海地下水模擬的準(zhǔn)確性. 由越流含水層地下水控制方程式(1)出發(fā)， 將潮汐波動(dòng)以正弦函數(shù)概化， 當(dāng)含水層系統(tǒng)距離海岸足夠遠(yuǎn)時(shí)， 地下水將逐漸趨于穩(wěn)定， 直至與平均海平面一致， 整個(gè)含水層系統(tǒng)在時(shí)間t=0時(shí)刻的初始地下水位與平均海平面相同， 當(dāng)t>0時(shí)， 自由含水層的地下水位依然保持恒定. 故潮汐邊界條件、 內(nèi)陸邊界條件和初始邊界條件表示如下

(26)

式中：A為潮汐波動(dòng)振幅；t0為潮汐波動(dòng)周期；Ln為越流含水層系統(tǒng)的離岸距離.

表1顯示各項(xiàng)計(jì)算參數(shù)的取值情況， 該案例空間邊界點(diǎn)n1離散為31點(diǎn)，每點(diǎn)間距100 m，時(shí)間邊界點(diǎn)n2離散為66點(diǎn)，每點(diǎn)間隔1 h，階數(shù)為25階. 沿海地區(qū)含水層地下水初水位水平的情況如圖5所示.

表1 沿海地下水模擬參數(shù)表

圖5 Trefftz法與有限差分法地下水位波動(dòng)曲線對比Fig.5 Comparison of groundwater level fluctuation curve between Trefftz method and finite difference method

圖5中直觀地反映了, 由Trefftz法計(jì)算得到的波動(dòng)曲線與后項(xiàng)有限差分法計(jì)算得到的結(jié)果基本一致， 能夠進(jìn)行密切擬合， 在進(jìn)一步證明Trefftz法于沿海地下水模擬研究適用的同時(shí)， 也直觀地描述了沿海地區(qū)地下水位的波動(dòng)特性. 由于受潮汐波動(dòng)的影響， 不僅與潮汐具有相似的正弦函數(shù)波動(dòng)特性， 始終在海平面附近上下波動(dòng)， 還與潮汐波動(dòng)具有相同的周期. 同時(shí)隨著潮汐波不斷向內(nèi)陸傳播， 近海水位變化比遠(yuǎn)岸水位變化更加明顯, 直至離岸距離為1.5 km， 遠(yuǎn)岸水位的波動(dòng)振幅已經(jīng)逐漸穩(wěn)定， 趨近于零.

3 結(jié)語

1) 針對越流含水層系統(tǒng)中承壓地下水建立模型， 無需使用傳統(tǒng)的時(shí)間推算法就可以解決暫態(tài)問題. 通過數(shù)值驗(yàn)證表明， Trefftz法的數(shù)值結(jié)果不僅與解析解高度吻合， 而且能夠產(chǎn)生極高精度的數(shù)值解， 從而證明該方法可便捷地應(yīng)用于解決越流含水層系統(tǒng)地下水問題.

2) 沿海地區(qū)地下水的演化規(guī)律對許多濱海工程有著重要影響， 以越流含水層系統(tǒng)的概念模型為基礎(chǔ)， 利用Trefftz法研究沿海地區(qū)越流含水層地下水的水位變化和波動(dòng)， 并且和有限差分法的計(jì)算結(jié)果進(jìn)行比較， 二者結(jié)果均反映了沿海地區(qū)地下水位隨時(shí)間的波動(dòng)和離岸距離的增大而逐漸衰減的高度擬合情況. 驗(yàn)證了Trefftz法應(yīng)用于沿海越流含水層的可行性與高效性. 因此， Trefftz法的適用范圍可能會(huì)在不久的將來擴(kuò)展到沿海工程問題.