時(shí)間分?jǐn)?shù)階的非飽和滲流數(shù)值分析及其應(yīng)用*

朱帥潤(rùn)，李紹紅，鐘彩尹，吳禮舟

（1.重慶交通大學(xué) 山區(qū)橋梁及隧道工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，重慶 400074；2.上海交通大學(xué) 土木工程系，上海 200240；3.西南交通大學(xué) 地質(zhì)工程系，成都 610031）

引 言

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

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

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

1 時(shí)間分?jǐn)?shù)階Richards 方程的數(shù)值近似

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

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

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

對(duì)于式（2）的右側(cè)分?jǐn)?shù)階導(dǎo)數(shù)項(xiàng)需采用Caputo 分?jǐn)?shù)階導(dǎo)數(shù)定義[21]：

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

式（6）可以進(jìn)一步簡(jiǎn)化為

對(duì)式（2）的左側(cè)進(jìn)行有限差分離散，有

其中

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖3顯示了t=2 h 時(shí)不同SWCC 形式下所計(jì)算的壓力水頭曲線，其中Gardner 和FX 模型的計(jì)算結(jié)果十分接近，在同一深度時(shí)兩者均小于VGM 模型的計(jì)算值.在同一深度可以發(fā)現(xiàn)當(dāng)γ＜1時(shí)，壓力水頭數(shù)值隨著γ的減小有減小的趨勢(shì)；當(dāng) γ＞1時(shí)，壓力水頭數(shù)值則隨著 γ的增大有增大的趨勢(shì).此外，不同形式下階次 γ對(duì)壓力水頭曲線的影響也存在明顯差異，在Gardner 和FX 模型中壓力水頭值隨著 γ的變化有較大的增幅（圖3（a）、（b）），而在VGM 模型中壓力水頭值隨著γ的變化增幅較?。▓D3（c））.這個(gè)結(jié)果表明了分?jǐn)?shù)階階次γ 相對(duì)于1 越大，入滲速度越快，階次γ 相對(duì)于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 入滲實(shí)驗(yàn)與應(yīng)用

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

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

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

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

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

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

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

數(shù)值模擬中假設(shè)模型深度為60 cm，總模擬時(shí)間選擇為160 h，空間離散步長(zhǎng) Δz=0.5 cm，時(shí)間步長(zhǎng)為0.1 h.圖7是采用VGM 和FX 模型得到的瞬態(tài)孔隙水壓力的計(jì)算結(jié)果，當(dāng)γ=1 時(shí)，兩種模型描述的經(jīng)典Richards 方程所獲得的數(shù)值解與實(shí)測(cè)值均存在較大偏差，其中FX 模型得到的數(shù)值解（圖7（c））相較于VGM 模型得到的數(shù)值解（圖7（a））偏差更大；當(dāng) γ=0.97 時(shí)，VGM 模型描述的時(shí)間分?jǐn)?shù)階Richards 方程所獲得的數(shù)值解可以與實(shí)測(cè)值吻合得很好（圖7（b）），如表3所示，其均方根誤差（RSE）和相對(duì)誤差（RE）均有很大程度的降低.t=24 h時(shí)，F(xiàn)X 模型描述的時(shí)間分?jǐn)?shù)階Richards 方程的RSE 和RE 均小于 γ=1的數(shù)值，也就是擬合效果有所提高，但其數(shù)值解與實(shí)測(cè)值仍然存在較大偏差（圖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í)間分?jǐn)?shù)階Richards 方程中階次γ的取值可以在區(qū)間[0.95,1.05]中選擇.結(jié)果表明，實(shí)驗(yàn)中出現(xiàn)的反常擴(kuò)散現(xiàn)象用VGM 模型的時(shí)間分?jǐn)?shù)階Richards 方程來(lái)進(jìn)行描述是更加合理的，與實(shí)測(cè)數(shù)據(jù)具有更好的擬合效果.

3 結(jié) 論

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

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

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