羅曉輝，朱帥潤，陳驕銳，吳禮舟,2

(1.成都理工大學(xué)環(huán)境與土木工程學(xué)院，地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，四川成都，610059；2.重慶交通大學(xué)省部共建山區(qū)橋梁及隧道工程國家重點(diǎn)實(shí)驗(yàn)室，重慶，400074)

全球變暖增加了諸如暴雨等極端天氣出現(xiàn)的頻次，導(dǎo)致了眾多的非飽和/飽和土淺層滑坡的地質(zhì)災(zāi)害[1?4]。因此，降雨入滲的模擬對邊坡穩(wěn)定性分析、滑坡地質(zhì)災(zāi)害等預(yù)測具有重要意義[5?10]。RICHARDS 方程[11]常被用以描述非飽和土降雨入滲問題，其中滲透系數(shù)與土?水特征曲線息息相關(guān)，獲得該方程的解析解較困難[12?14]，因此，地下水滲流問題常采用數(shù)值方法來近似求解[15?16]。關(guān)于降雨入滲問題的數(shù)值解或解析解[17?19]已有大量研究，例如：SRIVASTAVA等[20]使用指數(shù)函數(shù)形式表示的土?水特征曲線獲得了一維均質(zhì)和兩層土瞬態(tài)壓頭分布的解析解；TRACY[21]基于線性化的RICHARDS 方程獲得了RICHARDS 方程二維和三維的解析解表達(dá)式；KU等[22]使用Gardner模型線性化了RICHARDS 方程；孫長帥等[19]通過指數(shù)函數(shù)來描述土?水特征曲線、滲透系數(shù)函數(shù)，得到了關(guān)于一維非飽和土入滲過程的RICHARDS 方程的解析解。此外，用于求解可變飽和滲流問題的一些先進(jìn)的數(shù)值方法在某種條件下表現(xiàn)出較高的數(shù)值精度、計(jì)算效率和易于實(shí)施的優(yōu)點(diǎn)，如非正交網(wǎng)格的有限差分法[23]、無網(wǎng)格方法[24]、元胞自動機(jī)方法等[25]。例如，DOLEJ?í等[26]提出了自適應(yīng)高階時(shí)空不連續(xù)Galerkin方法來解決非飽和流的問題，以獲得更高的效率和精度。WU 等[27]使用Chebyshev譜方法求解一維邊坡的RICHARDS 方程，與有限差分法相比具有更高的精度。

雖然不同數(shù)值方法求解RICHARDS 方程具有不同的優(yōu)點(diǎn)，但研究一種改進(jìn)的迭代方法來求解由上述數(shù)值方法離散獲得的線性方程組，以提高數(shù)值收斂率依然很有意義。陳曦等[18]提出了短項(xiàng)混合欠松弛法求解RICHARDS 方程，并對其實(shí)用性和可靠性進(jìn)行了驗(yàn)證。此外，由于預(yù)處理可極大地加快迭代方法的收斂速度，故預(yù)處理技術(shù)也被廣泛地用于各種迭代方法中[28]。LIU 等[29]提出預(yù)處理的共軛梯度法減小了病態(tài)線性系統(tǒng)的條件數(shù)并獲得了較好的收斂效果。ZHU 等[17]提出了一種改進(jìn)的雙邊平衡Jacobi迭代方法，可有效地求解層狀非飽和土的嚴(yán)重病態(tài)線性方程組。WANG 等[30]提出了預(yù)處理的改進(jìn)的Hermitian 和skew-Hermitian分裂迭代方法，可有效地求解分?jǐn)?shù)階非線性Schr?dinger方程。

Gauss-Seidel 迭代法是求解線性方程組的一種經(jīng)典并且常用的迭代方法，但是該方法的收斂速度較慢[31]。SOR 迭代法的收斂速度比Gauss-Seidel迭代法的收斂速度高，但是其松弛系數(shù)的選擇是較困難的[31]。然而，根據(jù)迭代矩陣的最小譜半徑可以確定最優(yōu)松弛系數(shù)，結(jié)合矩陣分裂和預(yù)處理技術(shù)，SOR 迭代法的收斂率可以得到進(jìn)一步提高，因此，提出了2 種改進(jìn)方式的SOR 迭代，即紅黑排序SOR 迭代法(RB-SOR)和采用對角尺度化預(yù)處理子的預(yù)處理SOR迭代法(P-SOR)。

本文作者首先采用GS 和SOR 及2 種改進(jìn)的SOR 迭代算法求解穩(wěn)態(tài)和瞬態(tài)的線性化RICHARDS 方程，討論在不同網(wǎng)格長度下預(yù)處理技術(shù)對系數(shù)矩陣條件數(shù)的影響情況，進(jìn)而討論4種方法的收斂速度，以及在不同網(wǎng)格長度和時(shí)間步長下改進(jìn)的SOR 方法相對于傳統(tǒng)方法(GS 和SOR)的加速情況；同時(shí)，在不同網(wǎng)格長度和時(shí)間步長下，對采用不同迭代方法獲得的數(shù)值解與解析解之間的均方根誤差進(jìn)行比較，以檢驗(yàn)其數(shù)值精度和穩(wěn)定性。

1 控制方程和離散化

RICHARDS 方程通常被用來描述非飽和土體中地下水滲流問題，其三維壓力水頭(h)形式的RICHARDS方程(RE)可以寫成[11]

式中：Kx(h)，Ky(h)和Kz(h)分別為相對于x，y和z方向的滲透系數(shù)函數(shù)；h為壓力水頭；C(h)為容水度。由于土水特征曲線是非線性的，容水度可以描述為

式中：θ為含水率。求解RICHARDS方程需要3個(gè)特征函數(shù)，即非飽和土滲透系數(shù)函數(shù)、土水特征曲線以及含水率函數(shù)。將式(2)代入式(1)，對于一維非飽和土垂直入滲的RICHARDS方程可簡化表示為

式中：z為垂直坐標(biāo)。對于非飽和土體中的滲透系數(shù)和含水量函數(shù)，GARDNER提出如下指數(shù)形式[14]：

式中：θs和θr分別為飽和含水量和殘余含水量；Ks為土體的飽和滲透系數(shù)；αg為土體的孔隙大小分布參數(shù)。為了線性化RICHARDS 方程，通常先將非飽和土的滲透率進(jìn)行歸一化，則相對滲透系數(shù)Kr可表示為

進(jìn)而，引入一個(gè)新的參數(shù)h*，定義為

式中：χ為常數(shù)，χ=eαghd；hd為干燥土的壓力水頭。因此，對式(7)相對于z軸進(jìn)行求導(dǎo)，可得

從式(7)可以看出0 ≤h*<1，因此，

進(jìn)而對式(5)求導(dǎo)，可表達(dá)為

因此，將式(8)～(10)代入式(3)，可獲得線性RICHARDS方程：

式中：c=αg(θs-θr)/Ks。

采用有限差分方法(FDM)來離散式(11)，獲得1 個(gè)近似解。對于空間離散化，使用中心差分格式；對于時(shí)間離散，采用隱式格式。因此，1D 線性化RICHARDS方程的有限差分格式可表示為

式中：i為沿z軸方向的離散節(jié)點(diǎn)；Δz為離散的空間步長；Δt為離散的時(shí)間步長；n為時(shí)間節(jié)點(diǎn)?？紤]穩(wěn)態(tài)情況時(shí)，其有限差分格式可簡化如下：

式(12)和式(13)所構(gòu)成的線性方程組可寫成矩陣形式：

式中：A，h*和f分別為N×N階方陣、N×1 列向量和N×1 列向量。式(14)可采用相關(guān)迭代法進(jìn)行求解，一旦獲得數(shù)值解h*，實(shí)際的壓力水頭可以通過如下公式換算：

2 方法和改進(jìn)

對于經(jīng)典的線性方程組的矩陣分裂迭代法如Jacobi 迭代法、Gauss-Seidel 迭代法以及松弛型迭代法[31]，矩陣A可分解為

式中：D為A的對角矩陣；U為A的上三角矩陣；L為A的下三角矩陣。Gauss-Seidel 迭代方法的迭代格式可表達(dá)為

式中：迭代矩陣G=(D-L)-1U；b=(D-L)-1f；k為迭代次數(shù)。

為了進(jìn)一步提高Gauss-Seidel 迭代法(GS)的收斂速度，引入松弛系數(shù)ω，定義ω為0～2，則迭代解可表達(dá)如下：

式中：uk+1表示為式(17)的迭代結(jié)果；因此，根據(jù)式(17)～(18)超松弛迭代法(SOR)的迭代格式可表示為

在式(19)中，ω的選擇通常是困難的，沒有一個(gè)通用的選取規(guī)則。因此，本文研究中采取迭代矩陣的最小譜半徑來定義最佳ω，進(jìn)而再進(jìn)行迭代求解。式(19)的迭代矩陣為

不同的松弛系數(shù)ω可以從區(qū)間(0，2)中選取，其取值間隔為0.02，然后根據(jù)選取的ω計(jì)算對應(yīng)的迭代矩陣譜半徑，進(jìn)而搜索最小譜半徑位置，得到最佳松弛系數(shù)ω。將獲得的ω代入迭代公式進(jìn)行求解，在迭代求解過程中，當(dāng)?shù)鷿M足時(shí)，則迭代停止，否則繼續(xù)。

2.1 紅黑排序的SOR迭代法(RB-SOR)

對離散節(jié)點(diǎn)重新編號，一種經(jīng)典的排序方法(紅黑排序)可將節(jié)點(diǎn)分裂為2個(gè)向量[31]，奇數(shù)點(diǎn)(紅點(diǎn))為一組h＊R，偶數(shù)點(diǎn)(黑點(diǎn))為一組h＊B，進(jìn)而基于紅黑分裂，離散后的線性方程組可表達(dá)為

式中：DR和DB為對角矩陣，它們的元素個(gè)數(shù)分別為紅和黑點(diǎn)數(shù)量。C和C1分別為與待求的黑、紅點(diǎn)相關(guān)的交互矩陣。h＊R和h＊B分別對應(yīng)紅點(diǎn)和黑點(diǎn)的待求解向量；f1和f2分別為差分方程解算區(qū)域的邊值和方程的已知項(xiàng)。用SOR 迭代方法求解線性方程組(21)，有

表達(dá)為

式(23)是紅黑排序SOR 迭代一步的矩陣表達(dá)式，實(shí)際計(jì)算時(shí)也可不直接用這些矩陣公式，可結(jié)合式(18)更簡便地獲得紅黑排序SOR 的數(shù)值解。具體地，先采用GS計(jì)算在紅點(diǎn)的值，然后獲得SOR的值：

利用，再采用GS 計(jì)算在黑點(diǎn)的值，然后得到SOR值：

2.2 預(yù)處理SOR迭代法(P-SOR)

雖然SOR迭代法相比于GS迭代有了很大的改善，但SOR 的收斂速度仍需進(jìn)一步提高。預(yù)處理技術(shù)被廣泛地用于處理病態(tài)線性方程組問題并獲得了較好的效果，尋求預(yù)處理子(M)應(yīng)具有的性質(zhì)[28]有如下幾點(diǎn)：

1)M在某種意義下是矩陣A的較好近似；

2)M的構(gòu)造所占用計(jì)算時(shí)間是非常小的，甚至可以被忽略；

3)預(yù)處理后的線性方程組比原方程組(式(14))的求解更容易。

因此，在研究中采用對角尺度化預(yù)處理子，這種預(yù)處理技術(shù)易于實(shí)施和計(jì)算，對稱超松弛(SSOR)預(yù)處理子定義為

當(dāng)用上述矩陣作為預(yù)處理子時(shí)，沒有必要像迭代法那樣去細(xì)致地選擇參數(shù)ω。可取ω=1 即得對稱Gauss-Seidel(SGS)預(yù)處理子，有

采用式(27)的預(yù)處理子對線性方程組(14)進(jìn)行左預(yù)處理，有

對式(28)進(jìn)行SOR迭代求解，獲得相應(yīng)的數(shù)值解，其迭代格式如下：

式中：DM，LM和UM為預(yù)處理后的矩陣；fM為預(yù)處理后的列向量。對于一個(gè)線性方程組Ax=b而言，可以采用條件數(shù)來評價(jià)給出的非奇異矩陣A是否病態(tài)，條件數(shù)定義如下[15]：

式中：矩陣范數(shù)‖ ‖為Frobenius范數(shù)。

綜上，根據(jù)RB-SOR 和P-SOR 方法原理編寫程序(MATLAB2014a)并計(jì)算求解，計(jì)算流程圖如圖1所示。

圖1 RB-SOR和P-SOR迭代方法的計(jì)算流程圖Fig.1 Calculating flow chart of RB-SOR and P-SOR iterative methods

為了驗(yàn)證提出方法的數(shù)值精度，采用均方根誤差(Rh)進(jìn)行比較，有

式中：和h*分別為分析解和數(shù)值解矩陣；j為節(jié)點(diǎn)數(shù)量；t為模擬時(shí)刻節(jié)點(diǎn)。為了更直觀地比較P-SOR 的加速效果，加速比定義SGS/P-SOR為(以GS為例)：

式中：TGS為采用GS 迭代法計(jì)算的CPU 時(shí)間；TP-SOR為采用P-SOR迭代法計(jì)算的CPU時(shí)間，其計(jì)算機(jī)的CPU為Intel i5-4590。

3 方法比較與驗(yàn)證

3.1 均質(zhì)非飽和土的一維穩(wěn)態(tài)數(shù)值解

第1個(gè)算例描述了均勻非飽和土中的一維穩(wěn)態(tài)地下水滲流，數(shù)學(xué)模型如圖2所示。

圖2 一維降雨入滲模型Fig.2 One-dimensional model of rainfall infiltration

假定土體厚度L為10 m，則飽和含水量θs、殘余含水量θr、土體的飽和滲透系數(shù)Ks和土的孔隙分布參數(shù)αg分別為0.35，0.14，1×10?4m/h 和8×10?3m?1[15]。上、下邊界條件分別可描述為：

令干燥非飽和土的壓力水頭為hd=-1000 m，根據(jù)式(7)，歸一化的邊界條件可寫為：

對于一維均質(zhì)非飽和土的穩(wěn)態(tài)流問題，TRACY[21]提出穩(wěn)態(tài)方程的解析表達(dá)式：

這個(gè)例子中收斂準(zhǔn)則ε設(shè)置為10?8。

在不同的Δz下，分別采用GS，SOR、RB-SOR和P-SOR迭代法進(jìn)行計(jì)算，結(jié)果如圖3和表1所示，圖3中：條件數(shù)是由矩陣A的范數(shù)乘積所表示，量綱為一，其數(shù)值越大，線性方程組的病態(tài)性越嚴(yán)重。由圖3和表1可知：通過預(yù)處理技術(shù)處理后，原系數(shù)矩陣A的條件數(shù)大幅減小，至少減小了1個(gè)數(shù)量級。P-SOR方法迭代次數(shù)是最小的，GS迭代法的迭代次數(shù)最大，在數(shù)值上后者約為前者的100倍。

圖3 不同網(wǎng)格大小下條件數(shù)的比較Fig.3 Comparison of condition number for different grid sizes

表1 求解穩(wěn)態(tài)RE的數(shù)值結(jié)果Table 1 Numerical results for solving steady-state RE

圖4所示為求解一維穩(wěn)態(tài)方程不同方法的收斂率的比較結(jié)果(Δz=0.10 m)。由圖4可見：GS 的收斂率最??；RB-SOR 的收斂率相較于SOR 有所提高，但提高幅度較小；P-SOR 的收斂率最快，并且相較于SOR有了更大提高。通過與解析解比較，上述4種迭代方法均有較高的精度，相對于解析解的絕對誤差，GS 能達(dá)到10?6的數(shù)量級，SOR 達(dá)到10?8的數(shù)量級，P-SOR 和RB-SOR 均達(dá)到10?9的數(shù)量級。圖5所示為一維穩(wěn)態(tài)地下水流的數(shù)值解和解析解的比較。由圖5可知：通過P-SOR法獲得的壓力水頭與解析解十分吻合，當(dāng)Δz分別為0.20，0.10和0.05 m時(shí)，其均方根誤差分別為1.478×10?9，1.439×10?9和1.277×10?9，這說明改進(jìn)的P-SOR 具有較高的數(shù)值精度和穩(wěn)定性。

圖4 求解一維穩(wěn)態(tài)方程不同方法的收斂率的比較(Δz=0.10 m)Fig.4 Comparison of convergence rates of different methods for solving one-dimensional steady-state equation(Δz=0.10 m)

圖5 一維穩(wěn)態(tài)地下水流的數(shù)值解和解析解的比較Fig.5 Comparison of numerical and analytical solution for one-dimensional steady-state groundwater flow

3.2 均質(zhì)非飽和土一維瞬態(tài)流數(shù)值解

第2個(gè)算例描述的是均質(zhì)非飽和土的一維瞬態(tài)地下水滲流，數(shù)學(xué)模型依然如圖2所示。除了非飽和土的孔隙大小分布參數(shù)αg=3.2×10-4m?1外，土體厚度L和其他模型參數(shù)與一維穩(wěn)態(tài)模型是一致的。邊界條件也與3.1節(jié)中描述的一致，初始條件為h(z,t=0)=-105m。對于均質(zhì)非飽和土中的一維瞬態(tài)地下水滲流問題，TRACY提出的瞬態(tài)流解析解表達(dá)式如下[21]：

式中：λm=kπ/L；μm=α2g/4+λ2m；t為模擬時(shí)間。

總模擬時(shí)間設(shè)置為10 h，迭代方法的收斂標(biāo)準(zhǔn)ε設(shè)置為10?8。為了驗(yàn)證所提方法的計(jì)算效率和精度，網(wǎng)格長度取0.02，0.04和0.05 m，時(shí)間步長取為0.0100，0.0050和0.002 5 h。

圖6所示為條件數(shù)隨不同的Δt和Δz的變化，由圖6可知：系數(shù)矩陣A與預(yù)處理后的系數(shù)矩陣M?1A的條件數(shù)隨網(wǎng)格長度的減小而變大，隨時(shí)間步長的減小而減小，然而，通過預(yù)處理方法，很好地降低了系數(shù)矩陣的條件數(shù)，可以看到M?1A的條件數(shù)均得到了大幅度減小，甚至接近1.0，這表明預(yù)處理過程能夠有效地降低矩陣的病態(tài)性。

圖6 條件數(shù)隨不同的Δt和Δz的變化Fig.6 Change of condition number with different Δz and Δt

表2所示為4種不同迭代方法在不同的網(wǎng)格長度和時(shí)間步長條件下的迭代次數(shù)和計(jì)算時(shí)間。由表2可知：在不同網(wǎng)格長度和時(shí)間步長下，相比GS 和SOR，RB-SOR 和P-SOR 方法均能降低迭代次數(shù)和減少計(jì)算時(shí)間，而P-SOR的性能比RB-SOR的更好。圖7所示為求解RE 不同方法的迭代次數(shù)和收斂率的比較結(jié)果(Δz=0.04 m)。由圖7可知：在不同的Δt下，P-SOR表現(xiàn)出最少的迭代次數(shù)，RBSOR 迭代次數(shù)比GS 和SOR 的少；GS 的收斂率最慢，RB-SOR 和P-SOR 的收斂率相較于SOR 都有所提高，而P-SOR的收斂率最高。

表2 求解瞬態(tài)流問題的數(shù)值結(jié)果Table 2 Numerical results for solving transient flow problem

圖7 求解RE不同方法的迭代次數(shù)和收斂率的比較(Δz=0.04 m)Fig.7 Comparison of iterations and convergence rates using different methods for solving RE(Δz=0.04 m)

圖8所示為通過P-SOR獲得的數(shù)值解和解析解的比較結(jié)果。由圖8可知：由P-SOR方法獲得的瞬態(tài)壓力水頭曲線與解析解十分吻合。圖9所示為t=10 h時(shí)不同迭代方法在不同網(wǎng)格尺寸下的Rh。由圖9可知：當(dāng)t=10 h時(shí)，Rh達(dá)到了10?5數(shù)量級，同時(shí)，不同迭代方法的均方根誤差隨著網(wǎng)格尺寸的減小而減小。相較于SOR，P-SOR和RB-SOR的精度隨著時(shí)間步長的減小都有小量提高。由于GS迭代法的均方根誤差較改進(jìn)的SOR迭代法大0.5～1個(gè)數(shù)量級，因此，此處并未納入比較范圍。結(jié)果表明，改進(jìn)的SOR 方法均有較強(qiáng)的數(shù)值穩(wěn)定性以及較高的計(jì)算精度(近似等價(jià)于SOR的精度)。

圖8 通過P-SOR獲得的數(shù)值解和解析解的比較Fig.8 Comparison of numerical solutions obtained by PSOR with analytical solutions

圖9 t=10 h時(shí)不同迭代方法在不同網(wǎng)格的RhFig.9 Rh of different iterative methods at different grid length when t=10 h

圖10所示為當(dāng)Δt=0.005 h 時(shí)，P-SOR 和RBSOR 相對于GS 和SOR 的加速比，由圖10可見：在不同的網(wǎng)格長度下，RB-SOR與P-SOR迭代法相對于GS 的加速比是最大的，SOR 的次之；當(dāng)Δt=0.005 h 時(shí)，2 種改進(jìn)方法相對于GS 的最小加速比分別為3.602 和3.891，相對于SOR 的最小加速比為1.914 和2.067。從數(shù)值結(jié)果來看，與于GS 和SOR相比，P-SOR與RB-SOR迭代法均表現(xiàn)出較高的計(jì)算效率，但是P-SOR 的加速效果比RB-SOR的加速效果好。

圖10 比較不同改進(jìn)的SOR相對于GS和SOR的加速比(Δt=0.005 h)Fig.10 Comparison sof speed-up ratios of different improved SOR relative to GS and SOR(Δt=0.005 h)

4 結(jié)論

1)在求解非飽和土一維降雨入滲問題時(shí)，與GS 和SOR 相比，2 種改進(jìn)的SOR 方法均可提高進(jìn)收斂率和提高數(shù)值精度。但P-SOR 比RB-SOR 具有更高的收斂率、更高的計(jì)算效率，這說明預(yù)處理技術(shù)相比較于矩陣分裂技術(shù)(紅黑排序)具有更好的性能，可更大程度地提高收斂速度。

2)在離散RICHARDS 方程時(shí)，形成的系數(shù)矩陣A的條件數(shù)往往并不理想，尤其在穩(wěn)態(tài)情況下條件數(shù)很大，因此，采用預(yù)處理技術(shù)有效減小系數(shù)矩陣A的條件數(shù)。在穩(wěn)態(tài)情況下，系數(shù)矩陣的條件數(shù)減少了至少1個(gè)數(shù)量級，在瞬態(tài)情況下，條件數(shù)能夠極大地減少，甚至接近于1.0。

3)基于這次研究對比，RB-SOR和P-SOR相比較于GS和SOR均有一定程度的加速效果，尤其在較細(xì)的網(wǎng)格條件下，2種改進(jìn)方法的加速效果十分顯著。相對于SOR，RB-SOR和P-SOR的加速比分別達(dá)43.93%和51.62%。因此，改進(jìn)的SOR尤其是P-SOR在地下水滲流模擬中有廣闊的應(yīng)用前景。