非飽和滲流模擬中非均勻空間網格的改進方法

朱帥潤，何 博，吳禮舟，李紹紅，卿毅偉

（1.上海交通大學土木工程系，上海 200240；2.重慶交通大學山區橋梁及隧道工程國家重點實驗室，重慶 400074；3.成都理工大學環境與土木工程學院，四川 成都 610059）

非飽和滲流問題廣泛存在于眾多領域，如地質工程中降雨作用下土質邊坡的穩定性分析[1?4]，地下水科學中溶質遷移模擬[5?6]，以及礦業工程中煤層注水和煤層氣開采等[7?8]。其中，Richards方程[9]是非飽和滲流數值模擬的基本方程，有效并可靠地數值求解Richards方程對相關領域的科學研究和生產具有十分重要的意義。

Richards方程通常有3種形式，包括壓力水頭形式、含水率形式和混合形式[10?12]。其中，壓力水頭形式的Richards方程備受關注。Richards方程的解析解可在一些簡化條件下獲得[13?14]。如，Parlange 等[15]求解推導了一維Richards方程的復雜級數解以便獲得入滲通量；Warrick等[16]假設了一個簡單的土水特征曲線并推導出一維Richards方程的解析解。然而，水力傳導系數和土水特征曲線在現實中非常復雜，這導致高度非線性的偏微分方程，因此解析解的獲得非常困難。因此，在實際條件下Richards方程的求解只能通過數值方法獲得[17?20]。

在Richards方程的數值求解中，數值離散方法通常是必要的，常用的空間數值離散方法包括有限差分法[21]（FDM）、有限體積法[22]（FVM）以及有限元法[23]。對于時間離散化，通常采用向后差分法[24]。例如，吳夢喜[25]發展了求解Richards方程具有更好數值穩定性的一般有限元算法。Zambra等[26]構造了在空間和時間上具有很高精確度的有限體積法，用于求解非線性Richards方程。Chávez-Negrete等[24]提出了改進的有限差分法并結合自適應步長的Crank-Nicolson方法用于求解Richards方程。雖然在數值方法上有很多改進，但傳統的FDM和FVM都需要非常精細的網格來解析復雜的滲流區域，從而導致更多的計算時間。因此，越來越多的研究關注于非結構網格的改善[27?28]。最近，Deng等[27]利用非正交網格的坐標變換便捷地求解了復雜地形中的飽和-非飽和滲流問題。Chávez-Negrete等[24]采用不規則的三角形網格加密邊坡表層，獲得了更可靠的數值解。Dolej?í等[28]開發了一種各向異性網格自適應技術，這種方法在不損失精度的情況下顯著減少了自由度的數量以及計算時間。鄭川東等[29]提出一種動態局部自適應空間網格方法用于求解水流水質耦合方程，相對于傳統結構均勻網格模型，該方法精度高、效率高，能夠自動捕捉計算敏感區域，準確模擬復雜流態和物質的遷移擴散。何金輝等[30]在離散元軟件PFC2D基礎上，開發了考慮流體動網格的顆粒-流體耦合算法，有助于提高模擬大變形下的一維固結試驗等案例的數值精度。此外，一些模擬軟件也采用了動態網格劃分方法以獲得更好的數值結果，如 ANSYS、FLUENT 等[31]。

由此可知，非均勻空間網格的使用在一定程度上可以有效降低計算量，同時也能保證計算精度。但在一些不利數值條件下，比如入滲于初始干燥土壤中，以及分層土上下兩層的水力傳導系數相差較大時，上述研究的非結構空間網格以及動態網格方法的實現往往較為復雜且不太合適，因此，本文提出采用一種非均勻Chebyshev空間網格改進的FDM數值離散過程，并與傳統均勻空間網格獲得的數值結果和解析解對比驗證。

1 數值求解過程

1.1 控制方程

通常，Richards方程可用于描述多孔介質中非飽和滲流問題，其一維壓力水頭形式可以表示為[32]：

式中：h——壓力水頭/m；

t——時間/h；

Kz(h)——相對于z方向的水力傳導系數/（m·s?1）；

C(h)——容水度/m?1，C(h)= ?θ/?h；

θ——含水率/%。

在非飽和土滲流問題中，本文假設水力傳導系數和含水率相對于壓力水頭的數學關系由下面的指數模型給出[33]：

式中：θr——殘余含水率/%；

θs——飽和含水率/%；

Ks——土體的飽和水力傳導系數/（m·s?1）；

α——土性相關的模型擬合參數。

1.2 均勻網格法

為了獲得式（1）的數值解，傳統的有限差分法采用均勻網格法進行數值離散[22]。首先，z軸上的間隔被分為N等分，模擬時間被分為M等分。進而，式（1）的一、二階空間導數離散采用中心差分格式，時間導數則采用后向差分格式，其有限差分格式為：

式中：i——沿z軸的離散節點（邊界節點除外）；

m——時間節點數；

?z——離散的均勻空間步長（圖1）；

圖1 均勻網格和Chebyshev網格的示意圖Fig.1 Schematic diagram of uniform grid and Chebyshev grid.

?t——離散的時間步長。

考慮穩態情況時，其有限差分格式可簡化如下：

式中：Ki+1/2、Ki?1/2——相鄰節點對應的水力傳導系數的調和平均值，可表示為：

式（4）（5）所構成的線性方程組可進一步寫成矩陣形式：

式中：A——(N?1)×(N?1)階三對角矩陣；

h——(N?1)維列向量；

b——(N?1)維列向量，此處向量b的首尾元素已包含了邊界條件。

1.3 改進的Chebyshev網格法

由于非飽和滲流問題中壓力水頭在邊界和土層分界面往往存在較大的變化，因此通常需要一個更精細的網格來進行加密。然而傳統的均勻網格在加密過程中，將產生過多的計算網格節點，計算費時甚至仍不夠精確?？紤]上述情況，提出一種Chebyshev網格坐標[34]：

式中：L——土層厚度。

如圖1（b）所示，Chebyshev網格節點只會在界面處進行高度加密，很大程度上降低了網格節點數量。結合Chebyshev網格的FDM離散格式可以表達為：

式中：?zi——Chebyshev網格節點之間的不等間距；

δzi——計算節點處前后Ki+1/2和Ki?1/2之間的不等間距，如圖1所示。

式（10）也可以簡化為如式（8）的矩陣形式，由于水力傳導系數和含水率的非線性關系，通常數值離散后需要采用非線性迭代法對系數矩陣A反復評估并迭代計算。其中，Picard法是比較經典和實用的非線性迭代方法，即：

式中：k——迭代次數。

對于迭代過程，當相對殘差滿足如下迭代容差時，迭代終止：

式中：ε——迭代容差，在這次研究中設置為10?8[35]。

基于上述Richards方程的Chebyshev網格離散格式，使用MATLAB（版本：R2014a）語言開發了有關非飽和滲流的程序。

2 數值評價

2.1 均質土的非飽和滲流

該案例描述的是均質非飽和土中的一維瞬態非飽和滲流[17]，假設土層厚度L= 10 m，指數模型參數設置為： α= 1×10?4、θs= 0.50、θr=0.11，以及飽和水力傳導系數Ks= 2.5×10?8m/s。此外，邊界條件可以寫為：

假設h0= ?105m，該瞬態非飽和流問題的精確解為[36]：

其中，

此外，選擇3個指標，即均方根誤差（RSE），相對誤差（RE）和最大相對誤差（MRE），驗證所提網格方法的計算精度：

式中：hi——模型數值計算值；

3個指標的值越小，表示模型數值計算的精度越高?？偰M時間設置為5 h，節點數分別取為100，150，200，時間步長（ ?t）分別取為 0.01，0.02，0.04 h。圖2（a）顯示了當節點數N=100時，不同時間步長下不同網格方法獲得的MRE，可以看出Chebyshev網格法獲得的MRE范圍在0.6%～4.5%之間，隨時間t的增加表現為先減小后增大，當t<4 h時，MRE隨著 ?t的減小而減??；而均勻網格法獲得MRE范圍在1.3%～49%之間，尤其當t>1 h時，MRE隨著t的增加而增大，并遠遠大于Chebyshev網格獲得的MRE。圖2（b）顯示了當時間步長 ?t=0.01 h時，不同離散網格節點數下不同網格方法獲得的MRE，可以發現均勻網格獲得MRE隨著時間t的增加而增大，當t>1 h時，其誤差越來越大；而Chebyshev網格法獲得的MRE范圍僅在0.3%～2.3%之間，隨著時間t的增加表現為開口向上的拋物線，此外，兩種方法的MRE隨著網格節點數N的增加均有減小的趨勢。

圖2 不同數值條件下不同網格方法的最大相對誤差的比較Fig.2 Comparison of the maximum relative error of different grid methods under different numerical conditions

如圖3所示，在 ?t=0.01 h和N=200條件下采用兩種網格方法獲得數值解，并與精確解進行比較?？梢园l現，均勻網格法獲得的數值解與精確解存在較大的偏差，尤其在時間t>2 h后，如圖3（a）所示；而Chebyshev網格法獲得數值解與精確解十分吻合，隨著時間t的增加，相對誤差較小，如圖圖3（b）所示。

圖3 不同網格方法下獲得的數值解和精確解的比較Fig.3 Comparison of the numerical solutions obtained by different grid methods and exact solutions

正如表1所示，在t=5 h時，可以看到兩種方法的RSE和均勻網格的RE均隨著N的增加而減小，而Chebyshev網格的RE隨著N的增加而增大，但從數值上可以看出Chebyshev網格的RSE與均勻網格法相差近100倍，而RE與均勻網格法相差10倍以上。

表1 t=5 h時的數值精度Table 1 Numerical accuracy at t=5 h

這個案例說明了提出的Chebyshev網格方法并不因網格節點數的限制而改善精度，也就是說，該方法可以在較少的節點數下獲得較高的數值精度，又具有較小的計算開銷。

2.2 分層土的非飽和滲流

利用分層土中非飽和穩態流對提出的改進網格方法進行進一步的驗證，其數學模型如圖4所示。兩層土的模型參數設置為θs= 0.35，θr=0.14， α = 8×10?3。此外，土層 1厚度（L1）和土層 2厚度（L2）均設置為5 m，邊界條件與上一個案例一致，其中初始條件h0=?103m。對于兩層土，Chebyshev網格節點坐標可以定義為：

圖4 兩層非飽和土的Chebyshev網格Fig.4 Chebyshev grid of two-layer unsaturated soil

式中：N1和N2——土層1和土層2所離散的節點數，本案例中N1和N2均假設為40。

在分層土中，由于不同土壤類型組合的影響，不同的組合對非飽和滲流存在很大影響。正如表2所示，不同土壤的飽和水力傳導系數是不同的。為了進一步驗證所提網格的適用性，假設土層1的飽和水力傳導系數（Ks1）為固定值 10?1m/s，土層 2逐漸從粗砂土轉變為細黏土，其水力傳導系數（Ks2）從 10?2m/s變化為 10?9m/s（表3）。

表2 飽和土壤的水力傳導系數Table 2 Hydraulic conductivity values for saturated soils

表3 工況1至8的水力傳導系數Table 3 Hydraulic conductivity for cases 1 to 8

從圖5（a）可以看出，均勻網格和Chebyshev網格獲得的數值解是近似的，均可以較好模擬出工況2的非飽和滲流規律。圖5（b）描述的是工況4的數值結果，土層1、土層2的飽和水力傳導系數的比值為104，均勻網格無法精細刻畫分界面處的壓力水頭，而Chebyshev網格獲得的數值解卻可以較精細地刻畫分界面處壓力水頭的變化。正如圖5（c）所示，Chebyshev網格在工況6時也能精細地刻畫出分界面處的壓力水頭變化規律。這個案例進一步說明了所提的Chebyshev網格方法在一些不利入滲條件下尤其在兩層土的飽和水力傳導系數相差甚大時，均能以較少的網格節點數獲得更可靠的數值解。

圖5 不同工況下的數值解比較Fig.5 Comparison of numerical solutions for different cases.

2.3 與軟件Hydrus-1D的對比

將提出的網格方法與軟件Hydrus-1D進行對比以進一步驗證改進網格方法的準確性。其中，土水特征曲線（SWCC）采用Van Genuchten模型[37]來描述：

式中：s——基質吸力/kPa，可近似為負孔隙水壓力；

Se——有效飽和度，Se=1/[1+(αvs)n]m；

αv、n、m——與土性相關的模型參數，并且m=1 ? 1/n。

采用甘肅省天水市的次生黃土進行土柱入滲實驗。通過現場入滲過程監測，飽和水力傳導系數為1.5×10?7m/s。如圖6所示，根據土壤含水率和基質吸力的實測數據，采用Van Genuchte模型擬合得到SWCC，其確定系數R2為0.95。此外，Van Genuchte模型的擬合參數如圖6所示。數值模擬中假設模型深度為60 cm，總模擬時間為160 h，初始條件設置為 ?260 kPa。在Hydrus-1D模擬中，空間離散采用均勻網格，其步長設置為0.5 cm?；诒疚牡木鶆蚓W格法和Chebyshev改進網格方法，時間步長和離散節點數分別設置為0.1 h和120。

圖6 Van Genuchten 模型擬合得到的SWCCFig.6 SWCC fitted by Van Genuchten model.

圖7顯示了采用Hydrus-1D軟件和本文提出方法的數值結果，可以看出提出方法的數值解與軟件Hydrus-1D的模擬結果比較吻合，而且Chebyshev網格獲得的數值解比均勻網格獲得的數值解更接近于Hydrus-1D的數值解。如表4所示，Chebyshev網格相對于Hydrus-1D數值解的RSE、RE和MRE均小于均勻網格的數值，這表明提出的改進網格方法在相同離散網格節點數下更加準確合理。

圖7 不同方法下獲得的數值解的比較Fig.7 Comparison of the numerical solutions obtained by different methods

表4 數值精度比較Table 4 Comparison of the numerical accuracy

3 結論

（1）非飽和滲流數值求解時，尤其在有限差分法中，數值離散條件的空間步長往往需要較小的步長才能獲得可靠的數值解。這往往使得計算費時，甚至在一些不利數值情況下數值精度也并不能得到很大改善。提出的Chebyshev網格方法在保持離散節點數不變的情況下，通過一個余弦函數在兩側界面進行加密，進而再數值求解獲得更可靠的數值解。

（2）數值結果表明，在極端干燥的初始條件下，所提Chebyshev網格方法獲得的數值解與精確解十分吻合；在兩層土中上下土層的水力傳導系數相差較大時，Chebyshev網格方法所獲數值解也能在較少的節點數下精細刻畫出壓力水頭在分界面處的變化規律，與常規的均勻網格方法相比，該方法可以在較少的節點數下獲得較高的數值精度，又具有較小的計算開銷。

（3）在分層土中由于界面處水力傳導系數的平均化，往往也會使得計算效率下降，本文提出的空間網格方法可以在降低離散節點數層面上提高計算效率。同時，對比軟件Hydrus-1D，可以發現提出的方法與Hydrus-1D的數值解是比較吻合的，并且相較于均勻網格方法，提出的Chebyshev網格方法在相同離散網格節點數下是更加準確合理的。