朱帥潤，李紹紅，鐘彩尹，吳禮舟

（1.重慶交通大學 山區(qū)橋梁及隧道工程國家重點實驗室，重慶 400074；2.上海交通大學 土木工程系，上海 200240；3.西南交通大學 地質工程系，成都 610031）

引 言

非飽和滲流問題廣泛存在于巖土工程、邊坡工程以及地下水科學等領域[1-5].經(jīng)典的Richards 方程[6]可以描述孔隙介質中的非飽和滲流過程.由于水力傳導系數(shù)與土水特征曲線的非線性特征，Richards 方程的解析解通常是較難獲得的[7-8]，因此，Richards 方程的數(shù)值求解得到了很大程度的發(fā)展[9-10].數(shù)值求解過程中，一般會對Richards 方程進行數(shù)值離散，通常采用的數(shù)值離散方法包括有限差分法[11-12]、有限元法[13]和有限體積法等[14].例如，吳夢喜[15]發(fā)展了求解Richards 方程具有更好數(shù)值穩(wěn)定性的一般有限元算法.Liu 等[11]采用有限差分法數(shù)值離散了分層非飽和土中的滲流方程.Zambra 等[16]構造了在空間和時間上具有很高精確度的有限體積法，用于求解非線性Richards 方程.通過上述數(shù)值方法獲得的線性方程組需要迭代求解，此時Picard 法是最為實用和簡便的[17-19].

近年來，很多學者從擴散現(xiàn)象的角度出發(fā)研究地下水在非飽和土體中的輸運過程，發(fā)現(xiàn)其并不滿足經(jīng)典的Fick 擴散定律，也就是反常擴散過程[20].反常非飽和滲流本質上也是一種非Markov 過程，數(shù)值計算時則需要考慮滲流過程中的時間相關性或空間相關性[21].最近，關于分數(shù)階方程已有很多研究和進展.例如，Pachepsky 等[22]提出了廣義的Richards 方程，相比于整數(shù)階方程，時間分數(shù)階Richards 方程可以有效地表明土壤水分運移現(xiàn)象中存在的記憶效應.Gerasimov 等[23]結合時間分數(shù)階擴散方程更好地擬合了燒制黏土磚吸水和在水泥砂漿中入滲的實驗數(shù)據(jù).Gerolymatou 等[24]對實驗得到的具有適當初始值和邊界值的分數(shù)階Richards 方程進行了數(shù)值求解，并將模型結果與硅土磚中的水分入滲實驗數(shù)據(jù)進行了對比，得到了較好的擬合結果.王睿等[25-26]采用conformable、Caputo 導數(shù)推導了含水率格式的分數(shù)階Richards 方程，并與文獻[27]的實驗數(shù)據(jù)進行了對比研究.上述研究大多是對含水率格式的Richards 方程求解并與磚塊水分入滲實驗數(shù)據(jù)對比，因此，結合土柱入滲實驗，分析非飽和入滲過程中孔隙水壓力隨時間的變化規(guī)律，能夠為相關的巖土、邊坡工程提供有力支撐.

本文結合Caputo 導數(shù)得到了具有更廣泛滲流意義的時間分數(shù)階Richards 方程，采用有限差分法進行數(shù)值離散并迭代求解，以及對分數(shù)階參數(shù)和土水特征曲線進行了敏感性分析.最后，結合土柱入滲實驗數(shù)據(jù)，并與經(jīng)典Richards 方程的數(shù)值結果對比，對提出的時間分數(shù)階Richards 方程的擬合效果進行了驗證.

1 時間分數(shù)階Richards 方程的數(shù)值近似

Richards 方程可以被用于描述多孔介質中的非飽和滲流問題，其沿垂直方向z的一維壓力水頭（h）格式的Richards 方程可以表示為[28-29]

式中，h為壓力水頭；K(h)為相對于z方向的水力傳導系數(shù)；C(h)為容水度，其定義為C(h)=?θ/?h.在考慮入滲過程中的時間相關性時，時間分數(shù)階Richards 方程的表達如下：

式中，γ為關于時間導數(shù)的階次.當γ=1時，式（2）即退化為標準形式下的Richards 方程（1）.首先，對式（2）的求解，我們采用有限差分法進行數(shù)值離散，做網(wǎng)格剖分，令η，Δz分別為時間和空間離散步長，關于模擬時間被分為M等分，關于z軸上的長度被分為N等分，即

對于式（2）的右側分數(shù)階導數(shù)項需采用Caputo 分數(shù)階導數(shù)定義[21]：

Γ(·)為Gamma 函數(shù).式（5）中出現(xiàn)的壓力水頭導數(shù)的積分可以直接用數(shù)值微分公式逼近，其推導如下：

式（6）可以進一步簡化為

對式（2）的左側進行有限差分離散，有

其中

聯(lián)立式（7）、（8）可得式（2）的離散格式如下：

式中，k為迭代步數(shù)，A為對稱的三對角矩陣，h和b均為列向量.對于Picard 法的迭代求解過程，首先需要在每次迭代時導出線性方程組（12），然后使用基本的求解方法（例如Gauss 消元法、共軛梯度法）求解線性方程組[30].求解線性方程組后，式（12）中的系數(shù)矩陣A需要使用新的解向量hm,k+1進行更新，進而再次求解新的線性方程組.最后，前后兩次解向量的相對誤差滿足如下迭代容差時，迭代過程終止：

ε為迭代容差，本文中均設置為10-8[31].

此外，有許多數(shù)學模型可以描述上述非飽和土中的水力傳導系數(shù)和含水率與壓力水頭之間的數(shù)學關系[32-35].其中，Gardner 提出了如下指數(shù)模型[32]：

式中，Ks為飽和水力傳導系數(shù)；θs和 θr分別為飽和含水率和殘余含水率；α為土性相關的模型擬合參數(shù).用van Genuchten[33]和Mualem[34]（VGM）模型來描述水力傳導系數(shù)和含水率是比較經(jīng)典的：

式中，有效飽和度Se=1/[1+(|αh|)n]m；α，n和m均為與土性相關的模型擬合參數(shù)，m=1-1/n.

Fredlund-Xing（FX）模型描述的土水特征曲線如下[35]：

2 數(shù)值分析與應用

2.1 參數(shù)敏感性分析

基于上述時間分數(shù)階Richards 方程的離散格式，使用MATLAB（版本：R2014a）開發(fā)了有關非飽和滲流的程序.為了驗證分數(shù)階階次對非飽和滲流的影響，考慮三種土水特征曲線（SWCC）用于模擬在水平均質非飽和土中降雨入滲的一維瞬態(tài)流.數(shù)學模型如圖1所示，土層的厚度假設為L=10 m，其中飽和的水力傳導系數(shù)設置為Ks=6.0 × 10-2m/h[36].模型的邊界條件如下：

圖1 均質非飽和土的一維入滲模型Fig.1 The 1D infiltration model for homogeneous unsaturated soil

式中，h0代 表初始干燥土壤的壓力水頭值，并假設除去邊界點的初始條件為h(z,t=0)=-10 m.

總模擬時間設置為2 h，空間步長 Δz=0.1 m，時間步長η=0.01 h，階次γ 分別取為0.5，0.8，1.0，1.2 和1.5.我們在此示例中選擇沙質土壤，采用了Lu 和Likos[37]的實驗數(shù)據(jù).使用三種數(shù)學模型對實驗數(shù)據(jù)進行了擬合，如圖2所示，可以看出三種模型擬合得到的SWCC 在趨勢走向上基本一致，Gardner 模型在高基質吸力值處擬合不太好，VGM 和FX的確定系數(shù)R2分別為0.98 和0.99，在高基質吸力值處FX 模型擬合得更好.擬合參數(shù)包括土型相關的模型參數(shù)α，n，m，飽和含水率和殘余含水率，如表1所示.

圖2 SWCC的擬合曲線Fig.2 Fitting curves of the SWCC

表1 擬合參數(shù)Table 1 Fitting parameters

圖3顯示了t=2 h 時不同SWCC 形式下所計算的壓力水頭曲線，其中Gardner 和FX 模型的計算結果十分接近，在同一深度時兩者均小于VGM 模型的計算值.在同一深度可以發(fā)現(xiàn)當γ＜1時，壓力水頭數(shù)值隨著γ的減小有減小的趨勢；當 γ＞1時，壓力水頭數(shù)值則隨著 γ的增大有增大的趨勢.此外，不同形式下階次 γ對壓力水頭曲線的影響也存在明顯差異，在Gardner 和FX 模型中壓力水頭值隨著 γ的變化有較大的增幅（圖3（a）、（b）），而在VGM 模型中壓力水頭值隨著γ的變化增幅較小（圖3（c））.這個結果表明了分數(shù)階階次γ 相對于1 越大，入滲速度越快，階次γ 相對于1 越小，入滲速度越慢.

圖3 不同模型和階次γ 下獲得的壓力水頭曲線： （a） Gardner 模型；（b） Fredlund-Xing 模型；（c） van Genuchten-Mualem 模型Fig.3 Pressure head curves obtained under different models and orders γ: （a） Gardner model; （b） Fredlund-Xing model; （c） van Genuchten-Mualem model

2.2 入滲實驗與應用

入滲實驗所用土壤為中國甘肅省天水市的次生黃土.土柱模型裝置和實物如圖4所示，降水模型箱為透明有機玻璃，有機玻璃柱高度為1.45 m，填土高度為1.25 m，外徑20 cm，壁厚8 cm，距離土表面9 cm，14 cm，19 cm，24 cm，29 cm，39 cm，49 cm，59 cm，79 cm的側壁分別開有9個小孔，用于埋設體積含水率和基質勢傳感器.在模型頂部距離設計填土高度頂部設計有一出水口，高度為2 cm，用于控制積水高度.供水裝置主要由供水箱和滴水控制閥組成.供水箱材料為有機玻璃，分為內(nèi)圓和外圓，高度均為10 cm，內(nèi)圓外徑20 cm，底部中心開有一個5 mm 圓孔，用于安裝滴水控制閥，外圓外徑30 cm，開有3 cm 排水管孔，試驗中通過外環(huán)來調(diào)節(jié)內(nèi)環(huán)水面，使供水箱水位始終保持同一高度，保證供水水壓一致.滴水控制閥為醫(yī)用輸液管，可通過輸液管控制閥調(diào)節(jié)恒定流量.

圖4 土柱模型及供水裝置Fig.4 The soil column model and the water supply device

水分測定裝置由美國DECAGON 公司生產(chǎn)的EC-5 型體積含水率傳感器（圖5（a））、MPS-6 型基質勢傳感器（圖5（b））以及EM50 數(shù)據(jù)采集儀（圖5（c））構成.填土時將傳感器探頭埋入設計深度的土體內(nèi)，并與數(shù)據(jù)采集儀進行連接，再通過USB 線連接到筆記本電腦，即可獲得含水率和基質勢數(shù)據(jù).為了評價所提時間分數(shù)階模型的擬合效果，選擇了兩個指標，即均方根誤差（RSE）和相對誤差（RE），分別表達為

圖5 水分測定裝置：（a） EC-5；（b） MPS-6；（c） EM50Fig.5 Moisture measuring devices: (a) EC-5; (b) MPS-6; (c) EM50

式中，hk為模型計算值，h*k為土柱實驗中實測值，nn為測量的節(jié)點數(shù).兩個指標的值越小，表示模型的擬合效果越好.根據(jù)實驗現(xiàn)場數(shù)據(jù)，模型的邊界條件可以近似表達為式（22）、（23），初始干燥土壤的孔隙水壓力h0=-260 kPa.通過現(xiàn)場入滲過程監(jiān)測，飽和的水力傳導系數(shù)為Ks=1.296 cm/d.根據(jù)實測含水率與基質吸力數(shù)據(jù)，采用VGM 和FX 模型擬合得到的土水特征曲線如圖6所示，其確定系數(shù)R2分別達到了0.95 和0.98，說明FX 模型具有更好的擬合效果.此外，兩種模型的相關擬合參數(shù)如表2所示.

圖6 不同模型擬合得到的土壤水分特征曲線Fig.6 Soil moisture characteristic curves fitted by different models

表2 不同模型的擬合參數(shù)Table 2 Fitting parameters of different models

數(shù)值模擬中假設模型深度為60 cm，總模擬時間選擇為160 h，空間離散步長 Δz=0.5 cm，時間步長為0.1 h.圖7是采用VGM 和FX 模型得到的瞬態(tài)孔隙水壓力的計算結果，當γ=1 時，兩種模型描述的經(jīng)典Richards 方程所獲得的數(shù)值解與實測值均存在較大偏差，其中FX 模型得到的數(shù)值解（圖7（c））相較于VGM 模型得到的數(shù)值解（圖7（a））偏差更大；當 γ=0.97 時，VGM 模型描述的時間分數(shù)階Richards 方程所獲得的數(shù)值解可以與實測值吻合得很好（圖7（b）），如表3所示，其均方根誤差（RSE）和相對誤差（RE）均有很大程度的降低.t=24 h時，F(xiàn)X 模型描述的時間分數(shù)階Richards 方程的RSE 和RE 均小于 γ=1的數(shù)值，也就是擬合效果有所提高，但其數(shù)值解與實測值仍然存在較大偏差（圖7（d））.

圖7 比較不同SWCC 下不同階次γ 獲得的數(shù)值解： （a） VGM 模型，γ=1；（b） VGM 模型，γ=0.97；（c） FX 模型，γ=1；（d） FX 模型，γ=0.97Fig.7 Comparison of the numerical solutions obtained by different orders γ under different SWCC: （a） VGM model,γ=1; （b） VGM model,γ=0.97;（c） FX model,γ=1;（d） FX model,γ=0.97

表3 t=24 h的數(shù)值精度Table 3 Numerical accuracy at t=24 h

另外，由于黃土的入滲規(guī)律通常接近于經(jīng)典Richards 方程，因此提出的時間分數(shù)階Richards 方程中階次γ的取值可以在區(qū)間[0.95,1.05]中選擇.結果表明，實驗中出現(xiàn)的反常擴散現(xiàn)象用VGM 模型的時間分數(shù)階Richards 方程來進行描述是更加合理的，與實測數(shù)據(jù)具有更好的擬合效果.

3 結 論

本文基于非飽和滲流的Richards 方程，從反常擴散角度引入時間分數(shù)階導數(shù)，結合Caputo 導數(shù)定義，得到了具有更廣泛滲流意義的時間分數(shù)階Richards 方程，并結合土柱入滲實驗數(shù)據(jù)與其數(shù)值解進行了對比，得到了如下結論：

1） 時間分數(shù)階Richards 方程具有廣泛的適用性，當分數(shù)階階次γ=1 時，方程退化為經(jīng)典的Richards 方程，可以描述經(jīng)典的滲流過程；當分數(shù)階階次不等于1 時，則可以描述具有反常擴散性質的滲流過程，其中，當γ＜1時，壓力水頭數(shù)值隨著 γ的減小而減小，非飽和滲流表現(xiàn)為次擴散，而當 γ＞1時，壓力水頭數(shù)值則隨著 γ的增加而增加，表現(xiàn)為超擴散.

2） 從數(shù)值結果來看，不同土水特征曲線形式對分數(shù)階階次敏感程度是不同的，其中Gardner 和FX 模型中所計算的壓力水頭曲線隨著階次 γ的變化有較大的增幅，而在VGM 模型中所計算的壓力水頭曲線隨著γ的變化有較小的增幅.通過土柱入滲實驗，進一步驗證了VGM 模型的時間分數(shù)階Richards 方程具有更好的擬合效果，能夠更好地描述地下水在非飽和土中的滲流過程.