蘇 燕, 洪黎丹, 張瓏騰, 顧承宇

(1.福州大學(xué)土木工程學(xué)院, 福建 福州 350108; 2.臺灣海洋大學(xué)河海工程系，臺灣 基隆 20224)

0 引言

作為世界上受地質(zhì)災(zāi)害影響最嚴(yán)重的國家之一，2016年我國共發(fā)生地質(zhì)災(zāi)害9 710起，其中滑坡7 403起，占地質(zhì)災(zāi)害總數(shù)的76.2%[1].前人研究表明，降雨入滲過程中，非飽和土邊坡負(fù)孔隙水壓力變化是導(dǎo)致其失穩(wěn)的最主要因素[2-3].孔隙水壓力不僅與降雨條件有關(guān)，也與坡面的土性參數(shù)、土層分布以及初始孔隙水壓力有關(guān).降雨入滲的過程將導(dǎo)致土體性質(zhì)干燥狀態(tài)的動(dòng)態(tài)差異，在數(shù)值計(jì)算時(shí)極易形成病態(tài)矩陣，使計(jì)算結(jié)果不易收斂.基于網(wǎng)格單元的傳統(tǒng)有限元法和有限差分法在模擬非飽和土滲流問題時(shí)，由于濕潤峰的動(dòng)態(tài)性，滲透系數(shù)以及飽和含水量的強(qiáng)非線性，容易使網(wǎng)格發(fā)生畸變和扭曲，常需要重構(gòu)網(wǎng)格來解決，導(dǎo)致計(jì)算效率降低，甚至出現(xiàn)不收斂的現(xiàn)象[4].

針對降雨入滲條件下的邊坡穩(wěn)定性問題，學(xué)者們采用不同方法模擬非飽和土壤中的降雨入滲規(guī)律.羅渝等[5]利用有限差分法研究非飽和層降雨滑坡的運(yùn)動(dòng)特征.榮冠等[6]采用有限元法分析非飽和層土壤的降雨入滲機(jī)理.張社榮[7]結(jié)合強(qiáng)度折減技術(shù)，對非飽和層邊坡的的瞬態(tài)滲流場進(jìn)行求解.土壤水分特征曲線是影響非飽和層滲流變化的主要因素，學(xué)者們利用不同的滲流模型分析非飽和層土壤的瞬態(tài)滲流過程.甘永德等[8]運(yùn)用Modified Green-Ampt模型，模擬土壤積水和非積水情況下的入滲過程.Cho等[9]基于Mein-Larson模型，提出了非飽和層濕潤鋒深度的計(jì)算方法.

徑向基函數(shù)進(jìn)行插值作為多變量逼近理論中的一種數(shù)值方法，具有形式簡單、Euclidean距離各向同性，數(shù)值模擬結(jié)果與空間維數(shù)無關(guān)等優(yōu)點(diǎn)[10].Hard等[11]提出了以多元二次曲面函數(shù)作為基函數(shù)的插值方法，即MQ插值法，用來求解不規(guī)則地形表面的二次函數(shù)方程.Kansa等[12]將離散數(shù)據(jù)擬合和函數(shù)近似的MQ插值法與配點(diǎn)型無網(wǎng)格法結(jié)合，用于處理拉普拉斯方程，稱為MQ徑向基函數(shù)配點(diǎn)法(multiquadric radial basis function collocation method，MQ RBFCM).這種方法現(xiàn)已在很多領(lǐng)域有著良好的應(yīng)用.蘇李君等[13]將MQ配點(diǎn)法與有限差分法相結(jié)合，用于數(shù)值模擬非飽和土壤的水分運(yùn)動(dòng).為進(jìn)一步研究邊坡降雨入滲的瞬態(tài)滲流問題，針對非飽和層土壤特性的強(qiáng)非線性，本研究利用MQ配點(diǎn)法，建立一種數(shù)值計(jì)算模型，通過與前人的解析結(jié)果進(jìn)行對比，驗(yàn)證其數(shù)值模型的準(zhǔn)確性與穩(wěn)定性.在此基礎(chǔ)上，研究降雨入滲對非飽和無限邊坡的滲流場影響.

1 數(shù)學(xué)模型

1.1 控制方程

Tracy[14]根據(jù)Gardner提出的土水特征曲線模型，將高度非線性的Richards方程進(jìn)行線性化處理，得到一維線性化Richards方程：

(1)

圖1 均質(zhì)無限長斜坡示意圖Fig.1 Homogeneous infinite slope schematic diagram

式中：h*=eαgh-χ為線性化后的壓力水頭；χ=eαghd，hd=10-5為非飽和土極度干燥時(shí)的壓力水頭；Kθ=ks/(θs-θr)，ks為飽和滲透系數(shù)，θr為殘余含水量，θs為飽和含水量；Kαθ=Kθ/αg，αg為土壤孔隙大小分布系數(shù).式(1)僅針對無限邊坡成立，若考慮邊坡傾角為γ的無限邊坡問題，定義如圖1所示的(x*-z*)坐標(biāo)系，兩坐標(biāo)之間的關(guān)系為：

(2)

將式(2)代入式(1)得到：

(3)

1.2 邊界條件

對于坡底邊界，其壓力水頭恒定為ψ0，坡面為定流量邊界，式中q為降雨入滲量即：

(4)

2 數(shù)值算法

2.1 MQ RBFCM

MQ RBFCM求解偏微分方程時(shí)，首先將區(qū)域離散成N個(gè)點(diǎn)，區(qū)域邊界點(diǎn)滿足邊界條件，內(nèi)部點(diǎn)滿足偏微分方程， 建立N個(gè)代數(shù)方程.然后將待求解方程以徑向基函數(shù)近似函數(shù)表示，將近似函數(shù)及各階導(dǎo)數(shù)直接代入邊界條件或偏微分方程，則可將偏微分方程轉(zhuǎn)化為代數(shù)方程組進(jìn)行求解N個(gè)待定系數(shù)[15].假定待求解偏微分方程問題為：

(5)

式中：L和B為微分算子.在控制域Ω內(nèi)布置m個(gè)節(jié)點(diǎn)P1,P2, …,Pm，邊界上?Ω布置n個(gè)節(jié)點(diǎn)Pm+1,Pm+2, …,PN，然后令：

(6)

(7)

令N個(gè)節(jié)點(diǎn)各自滿足式(7)，得到N條方程，聯(lián)立方程求解待定插值系數(shù)λj，然后代回式(6)得到方程u(x)的解.

2.2 控制方程及邊界條件的離散

圖2 布點(diǎn)示意圖Fig.2 Point diagram

沿z*方向在控制域內(nèi)布M個(gè)內(nèi)部點(diǎn)，上邊界下邊界各布1個(gè)邊界點(diǎn)，共布N=M+2個(gè)點(diǎn)，如圖2所示.利用配點(diǎn)，將壓力水頭函數(shù)h*(z*)以MQ RBF形式離散為：

(8)

將離散后的壓力水頭函數(shù)h*(z*)代入控制方程和邊界條件，得到離散后的控制方程和邊界條件為：

(9)

(10)

(11)

2.3 瞬態(tài)滲流問題的求解

對于瞬態(tài)滲流問題的初始條件，借鑒前人方法，對上邊界施加一個(gè)較小的降雨強(qiáng)度qa，對邊坡進(jìn)行穩(wěn)態(tài)分析來求得初始孔隙水壓.根據(jù)M個(gè)內(nèi)部點(diǎn)滿足初始狀態(tài)下的穩(wěn)態(tài)控制方程，上下邊界點(diǎn)滿足各自的邊界條件，聯(lián)立可得N條方程，寫為矩陣形式：

{[A1]M×N, [A2]1×N, [A3]1×N}T{λ1,λ2,λ3, …,λM,λ1+M,λN}T={d1,d2,d3, …,dM,f1+M,gN}T

(12)

式中:

(13)

(14)

式中：H為系統(tǒng)矩陣，

(15)

根據(jù)已知初始狀態(tài)和邊界條件，可通過下式得到一組滿足條件的待定插值系數(shù)λj：

(16)

(17)

(18)

(19)

3 數(shù)值結(jié)果與討論

3.1 均質(zhì)無限邊坡

1) 計(jì)算模型與參數(shù).采用文[17]的無限邊坡，如圖1所示.斜坡土層的垂直厚度為2 m，斜坡水平長度為40 m，傾角γ=30°，飽和含水率θs=0.45，殘余含水率θr=0.15，飽和滲透系數(shù)ks=1.0×10-4cm·s-1，土壤孔隙大小分布系數(shù)α=0.01 cm-1.

2) 邊界條件.前期雨強(qiáng)為qa=2.8×10-11cm·s-1.邊坡上邊界邊界條件為降雨定流量邊界，降雨強(qiáng)度為qb.下邊界為定水頭邊界，ψ0=-1 m.

3) 計(jì)算結(jié)果.當(dāng)降雨強(qiáng)度qb為6.0×10-5cm·s-1小于ks和qb為3.0×10-4cm·s-1大于ks時(shí)，MQ RBFCM所得結(jié)果與前人的數(shù)值模擬及解析解所得的結(jié)果，在不同時(shí)刻的孔隙水壓力沿斜坡面垂直方向的分布情況如圖3所示.MQ RBFCM所得結(jié)果與前人的數(shù)值模擬及解析解吻合良好，表明運(yùn)用MQ配點(diǎn)法求解無限邊坡瞬態(tài)滲流問題是切實(shí)可行的.當(dāng)降雨強(qiáng)度大于飽和滲透系數(shù)時(shí)，李寧等[17]利用解析解求得的積水時(shí)間為11.297 h，本研究計(jì)算求得的積水時(shí)間為11.260～11.265 h，兩者的計(jì)算結(jié)果是比較接近的.

3.2 成層無限邊坡

1) 計(jì)算模型與參數(shù).采用文[16]的無限邊坡, 如圖4所示.成層斜坡中層1和層2的垂直厚度分別為1 m，斜坡水平長度為40 m，傾角γ=30°，如圖5所示.計(jì)算中層1和層2土體采用的非飽和土參數(shù)分別為： 飽和滲透系數(shù)ks1=3×10-3cm·s-1，ks2=3×10-4cm·s-1, 土壤孔隙大小分布系數(shù)α1=α2=0.01 cm-1，飽和含水率均為0.45，殘余含水率均為0.15.

圖3 不同時(shí)刻斜坡內(nèi)孔隙水壓力分布Fig.3 Distribution of pore water pressure at slopes in different time

圖4 雙層斜坡示意圖Fig.4 Sketch of double slope

圖5 雙層斜坡內(nèi)孔隙水壓力分布Fig.5 Pore water pressure distribution in double-layer slope

2) 邊界條件.qa=3.0×10-3cm·s-1為獲得初始負(fù)孔隙水壓力分布的前期雨強(qiáng)，qb=1.08×10-3cm·s-1為后期降雨強(qiáng)度.模型上邊界為定流量邊界，下邊界為定水頭邊界h=-1 m.

3) 計(jì)算結(jié)果.詹良通等[16]曾推導(dǎo)出成層斜坡降雨入滲的解析解.本研究參照文[16]，利用MQ RBFCM對該問題進(jìn)行求解，其結(jié)果對比前人解析解和有限元數(shù)值模擬計(jì)算結(jié)果可以發(fā)現(xiàn)，基本吻合如圖5所示，說明了運(yùn)用MQ RBFCM求解雙層斜坡降雨入滲孔隙水壓變化的可行性和正確性.

不同土壤滲透系數(shù)差異較大，變化范圍為1～10-8m·s-1，將土層1改為滲透系數(shù)為0.1 cm·s-1的砂土，改變土層2土壤，其滲透系數(shù)為10-2～10-9cm·s-1.傳統(tǒng)數(shù)值方法在求解互層土滲透系數(shù)差異較大時(shí)容易出現(xiàn)病態(tài)題，通過MQ RBFCM方法其條件數(shù)控制在1012以下, 如圖6所示，說明本方法求解系統(tǒng)較為良態(tài).將層1滲透系數(shù)定為10-1cm·s-1，層2為10-9cm·s-1，時(shí)間間隔為15 d, 歷時(shí)90 d.孔壓如圖7所示, 與實(shí)際物理現(xiàn)象相符，說明MQ RBFCM方法能較有效解決互層土滲透系數(shù)的極端差異在傳統(tǒng)數(shù)值方法中的病態(tài)問題.

圖6 案例1至8的條件數(shù)Fig.6 Number of conditions for Cases 1 to 8

圖7 雙層斜坡內(nèi)孔隙水壓力分布Fig.7 Pore water pressure distribution in double-layer slope

4 結(jié)語

應(yīng)用MQ 徑向基函數(shù)配點(diǎn)法求解無限長斜坡在雨強(qiáng)小于及大于土體飽和滲透系數(shù)的情況，及滲透系數(shù)極端差異的互層土體的的滲流問題.該結(jié)果與解析解對比來看，本方法得到的孔隙水壓力與解析解得到的水頭相吻合.說明本方法對求解無限長斜坡滲流問題是正確的，可用于預(yù)測不同降雨強(qiáng)度下斜坡內(nèi)孔隙水壓力隨時(shí)間的變化規(guī)律，以及非飽和土邊坡穩(wěn)定性分析.