路基邊坡雨水浸透變形模擬研究

馮振華，王利超，韓玉龍

摘要：為了減少由于雨水滲透引起的路基邊坡局部變形和失穩(wěn)現(xiàn)象，在概述模型實(shí)驗(yàn)現(xiàn)象的基礎(chǔ)上，采用飽和及非飽和滲流分析方法以及彈塑性有限元分析方法，模擬路基邊坡模型中水位的變化和滑動(dòng)破壞發(fā)生前的初始變形，建立漸進(jìn)式邊坡失穩(wěn)的預(yù)測(cè)方法，通過實(shí)驗(yàn)進(jìn)行再現(xiàn)分析，以明確邊坡初始變形機(jī)理。研究結(jié)果表明：1）在飽和及非飽和滲流模擬中，當(dāng)滲流系數(shù)設(shè)置值比滲透實(shí)驗(yàn)結(jié)果大時(shí)，模型實(shí)驗(yàn)的再現(xiàn)性良好;2）在彈塑性有限元分析中，通過邊坡水位升高導(dǎo)致變形系數(shù)（剛度）減小，進(jìn)而能夠分析出模型實(shí)驗(yàn)中邊坡的初始變形。研究過程為分析降雨期間路基邊坡穩(wěn)定性提供了新方法，其結(jié)果可為工程施工過程中避免發(fā)生大規(guī)模崩塌現(xiàn)象采取適當(dāng)措施提供依據(jù)，為施工標(biāo)準(zhǔn)的設(shè)計(jì)和制定提供重要參考。

關(guān)鍵詞：巖土工程;降雨模型實(shí)驗(yàn);初始變形;飽和及非飽和滲流分析;彈塑性有限元分析

中圖分類號(hào)：TU45文獻(xiàn)標(biāo)識(shí)碼：ADOI： 10.7535/hbgykj.2021yx05004

Simulation research on rainwater infiltration deformation of

roadbed slope

FENG Zhenhua，WANG Lichao，HAN Yulong

（Henan Transportation Planning and Design Institute Company Limited，Zhengzhou，Henan 450000，China）

Abstract：In order to reduce the local deformation and instability of subgrade slopes caused by rainwater infiltration，on the basis of an overview of the model experiment phenomena，the saturated and unsaturated seepage analysis and the elastoplastic finite element analysis were used to simulate the change of water level in the roadbed slope model and the initial deformation before sliding failure.Then the progressive slope instability prediction method was established.Finally，through the reconstruction analysis of the rules with experiments，the initial deformation mechanism of the slope was clarified.The results show that：1） in saturated and unsaturated seepage simulations，when the seepage coefficient is set larger than the result of the seepage test，the reproducibility of the model experiment is good;2） in the elastoplastic finite element analysis，the deformation coefficient （stiffness） can be reduced by increasing the water level of the slope，thereby the initial deformation of the slope in the model experiment can be analyzed.The study provides a new way for analyzing the stability of roadbed slopes during rainfall，and the results can provide a reference for taking appropriate construction measures to prevent the occurrence of large-scale collapse，which is also beneficial for formulating future construction standards.

Keywords：geotechnical engineering;rainfall model experiment;initial deformation;aturated and unsaturated seepage analysis;lastoplastic finite element analysis

河南地區(qū)每年7—9月為雨季，近年來大部分地區(qū)中雨、大雨以及暴雨日數(shù)呈增加趨勢(shì)，自然災(zāi)害頻繁。極端降水會(huì)引起路基土體結(jié)構(gòu)變形、坍塌，對(duì)公路路基影響巨大，嚴(yán)重影響交通運(yùn)輸安全。

豫東地區(qū)土質(zhì)多為粉質(zhì)砂土，降雨過程中邊坡土體結(jié)構(gòu)破壞的狀況主要通過模型實(shí)驗(yàn)進(jìn)行分析，確定邊坡的穩(wěn)定性[1-3]。對(duì)降雨過程中邊坡水位變化和崩塌機(jī)理的研究發(fā)現(xiàn)，粉質(zhì)砂土邊坡的穩(wěn)定性分析大多基于完全剛塑性體的極限平衡假設(shè)。但在全尺寸模型[4]和離心機(jī)模型[5]實(shí)驗(yàn)中，使用由粉質(zhì)砂土模型邊坡進(jìn)行降雨滑坡實(shí)驗(yàn)，隨著降雨造成的邊坡中的水位和孔隙水壓力的增加，邊坡會(huì)在某時(shí)間點(diǎn)發(fā)生局部變形，之后破壞逐漸向邊坡上部延伸，最終導(dǎo)致土體大規(guī)模滑動(dòng)，即漸進(jìn)式邊坡失穩(wěn)。此類研究中的邊坡破壞形態(tài)并非極限平衡法中所假設(shè)的瞬間滑動(dòng)。在漸進(jìn)式破壞發(fā)生前，坡面還出現(xiàn)了裂紋和局部變形。考慮到實(shí)際邊坡狀況，需要采取適當(dāng)?shù)氖┕し桨敢员阍诔跗谧冃吻胺乐雇馏w大規(guī)模滑動(dòng)崩塌發(fā)生。因此，闡明初始變形的機(jī)理，建立邊坡漸進(jìn)式破壞的預(yù)測(cè)方法顯得尤為重要。

近年來SPH法等粒子分析方法在對(duì)邊坡漸進(jìn)式失穩(wěn)的分析中被證明是有效的，并得到了坍塌滑動(dòng)過程中的影響因素分布的結(jié)論[6-8]。但在粒子分析方法中，需要的參數(shù)更多并且難以從現(xiàn)場(chǎng)實(shí)驗(yàn)勘測(cè)和室內(nèi)土力學(xué)實(shí)驗(yàn)中獲得。因此，通過有限元法簡(jiǎn)單地再現(xiàn)邊坡漸進(jìn)式失穩(wěn)具有更大的實(shí)際應(yīng)用價(jià)值。

本研究對(duì)模型邊坡實(shí)施了人工降雨實(shí)驗(yàn)，記錄了降雨過程中坡體的飽和度和孔隙壓力變化，并進(jìn)一步通過飽和及非飽和滲流模擬再現(xiàn)了模型邊坡中的水位變化，并嘗試通過彈塑性有限元法再現(xiàn)粉質(zhì)砂土邊坡模型中水位上升引起的邊坡初始變形來研究其整體變形機(jī)理。在彈塑性有限元分析中，主要針對(duì)巖土材料的變形系數(shù)展開研究，并假設(shè)邊坡地下水位升高時(shí)，變形系數(shù)隨有效應(yīng)力的減小和應(yīng)變的增加而產(chǎn)生折減。

1降雨邊坡模型實(shí)驗(yàn)

1.1模型概要

邊坡模型的概略圖如圖1所示。

水平支撐地層厚為70 mm，在其上制作了坡度比為1∶1.5，高600 mm、寬600 mm的邊坡模型，并在土箱壁上涂抹硅脂，從而減小摩擦。模型中使用的土體材料取自河南省道323高家溝大橋附近的粉質(zhì)砂土（土粒相對(duì)密度s=2.683 g/cm3，平均粒徑D50=0.165 mm，不均勻系數(shù)Uc=2.59，曲率系數(shù)U'c=1.32，細(xì)粒含量（質(zhì)量分?jǐn)?shù)）Fc=10.8%，最大干密度dmax=1.533 g/cm3，最優(yōu)含水率w = 19.5%），制作方法為普羅克特壓實(shí)法，干密度d /最大干密度dmax×100%=95%，w=19.5%。以邊坡體壓實(shí)度Dc=75%，w=13.0%為目標(biāo)，每次鋪設(shè)50 mm厚度土層，然后用振動(dòng)器壓實(shí)后再進(jìn)行下一層土層的鋪設(shè)。研究中，在坡箱底部放置一層平均直徑約10 mm的礫石材料，以防止由于水位快速上升發(fā)生滲透崩塌。

在邊坡模型的坡頂安裝位移計(jì)，在坡體中安裝孔隙水壓計(jì)和土壤水分儀用于測(cè)量土體含水量變化，儀器安裝布置如圖1所示。降雨量設(shè)置為40 mm/h，降雨持續(xù)450 min，其后采用灑水降雨。在支撐土層的3個(gè)位置安裝壓力計(jì)，以測(cè)量邊坡水位。

1.2數(shù)據(jù)分析

土體壓力計(jì)讀數(shù)隨降雨時(shí)間的變化如圖2所示，壓力計(jì)值表示了邊坡中地下水位高度。降雨持續(xù)450 min后直至模型邊坡崩塌過程中，隨著降雨不斷進(jìn)行，壓力計(jì)值不斷增大，說明邊坡模型內(nèi)土體水位持續(xù)上升;降雨551 min后，坡面土體開始發(fā)生變形。實(shí)驗(yàn)過程發(fā)現(xiàn)，壓力計(jì)值和孔隙水壓力值增加不同步（如圖3所示），這是因?yàn)閴毫τ?jì)值是在測(cè)點(diǎn)周圍的土體高度飽和后其讀數(shù)才開始上升，而孔隙水壓力在降雨過程中土體連續(xù)滲透作用下就會(huì)產(chǎn)生，并不是在土體高度飽和狀態(tài)下才會(huì)發(fā)生。由于壓力計(jì)讀數(shù)與孔隙壓力上升過程的趨勢(shì)基本一致，因此可判斷壓力計(jì)讀數(shù)是合理的。模型邊坡從輕微變形到坍塌的過程如圖4所示。圖中①處，在深度約200 mm的地方坡體發(fā)生了初始裂縫和輕微變形;②處，地下水開始從模型坡體產(chǎn)生裂縫的地方滲出;③處，一定時(shí)間后，滲水點(diǎn)上方的斜坡表面開始出現(xiàn)裂縫;④處，以坡體頂部裂縫為邊界的滑塌持續(xù)發(fā)生。另外從圖4中可以觀察到坡面下端的塌陷土變得泥濘。

自坡面有微小初期變形開始，邊坡截面形狀隨時(shí)間的變化如圖5所示。在初期變形發(fā)生10 min后邊坡下端發(fā)生了較大的變形，并且發(fā)生了小規(guī)模的崩塌。隨著崩壞規(guī)模逐漸變大，在坡面發(fā)生變形60 min后崩塌范圍延伸到了邊坡的頂部，最終形成大規(guī)模圓弧狀滑動(dòng)。

圖3中亦顯示了飽和度Sr和孔隙壓力增量Δuw隨著時(shí)間的變化，飽和度是根據(jù)土壤水分儀獲得的體積含水量以及與制作模型時(shí)的相同ρd（ρd=1.15 g/cm3）計(jì)算出。Δuw為澆灌過程中孔隙水壓的增量，其初始值為澆灌開始時(shí)的孔隙水壓值。

飽和度在開始灑水降雨的初期階段，邊坡底端（M14和M15）土體的Sr值有所上升，之后此處Sr值約為0.8，高于其他區(qū)域。這是因?yàn)榻?rùn)面在斜坡底部較快到達(dá)防滲層（在本研究中防滲層為支撐上部邊坡的堅(jiān)實(shí)土層），隨后形成水位使飽和度Sr值增加。從Δuw的變化可以看出，在降雨75 min左右，斜坡末端（P07）的Δuw上升，隨著降雨的滲透，Δuw趨向于從邊坡底端向模型中心（從P07到P01）方向上升。在模型底部飽和度Sr=0.8時(shí)，邊坡中間部分（M12）的Sr值增大，當(dāng)中部的Sr值和邊坡底端的Sr值大致相同時(shí)，發(fā)生初始變形。由此可知，造成邊坡初期變形的原因之一是由于邊坡底端和中部水位升高，而導(dǎo)致邊坡抗剪切強(qiáng)度下降[7]。

2灑水降雨時(shí)水位上升的模擬

2.1數(shù)值模型和條件

數(shù)值模型截面如圖6所示，數(shù)值模型邊界條件如圖7所示。

數(shù)值模型的尺寸與模型實(shí)驗(yàn)尺寸相同。滲流數(shù)值分析邊界條件：模型上部作為降雨入滲邊界，支撐層的左端作為排水邊界，模型的最底部和右緣均假定為不透水邊界。降雨條件設(shè)置與實(shí)驗(yàn)相同。數(shù)值分析的初始條件是當(dāng)安裝在邊坡和支撐層邊界之間的壓力值發(fā)生變化（降雨250 min），即支撐層的水位達(dá)到70 mm時(shí)，相當(dāng)于給支撐層底面施加70 mm的壓力水頭。與模型實(shí)驗(yàn)土體相同壓實(shí)度（Dc）試樣的浸透排水實(shí)驗(yàn)結(jié)果如圖8所示。

滲流分析中使用了水分特征曲線Van Genuchten模型[9]以及比滲透系數(shù)模型Mualen[10]。

θ=1+（α×ψ）n-m×（θs-θr）+θr，（1）

Kr=Se1/21-（1-Se1/m）m2，（2）

m=1-1/n，（3）

式中：θ為體積含水率;θr為殘留體積含水率;θs為飽和時(shí)體積含水率;ψ為基質(zhì)吸力，kPa;Se為有效飽和度;Kr為比滲透系數(shù)。

2.2數(shù)值分析

灑水降雨期間壓力計(jì)高度hm隨時(shí)間的變化如圖9所示。這里表示了使用與模型實(shí)驗(yàn)干密度相同的試樣進(jìn)行定水頭滲透實(shí)驗(yàn)得到的滲透系數(shù)k（k=klab）和使用與實(shí)驗(yàn)結(jié)果整合后逆向計(jì)算得到的滲透系數(shù)（k=kmod）的對(duì)比結(jié)果。

采用定水頭滲透實(shí)驗(yàn)結(jié)果的分析中，Ma01的數(shù)值結(jié)果始終大于實(shí)驗(yàn)結(jié)果，Ma02在t=400 min左右的再現(xiàn)性較高，之后數(shù)值結(jié)果遠(yuǎn)大于實(shí)驗(yàn)結(jié)果。在t=400 min之前，Ma03的數(shù)值結(jié)果略低于實(shí)驗(yàn)結(jié)果，此后數(shù)值解比Ma02中的實(shí)驗(yàn)結(jié)果大很多。是因?yàn)樗峁┑挠晁疀]有及時(shí)從模型邊坡中排出，從而導(dǎo)致模型中的水位被過高計(jì)算。通過使用表1所示kmod數(shù)值分析，可以很好再現(xiàn)實(shí)驗(yàn)結(jié)果（圖9中k=kmod）。

滲透坍塌部分的土體會(huì)從坡面流出?，F(xiàn)推測(cè)是由于土體壓實(shí)不均勻在坡體中產(chǎn)生了細(xì)小過水通道進(jìn)而導(dǎo)致局部排水，這些過水通道的出現(xiàn)，導(dǎo)致整個(gè)模型實(shí)驗(yàn)的滲透系數(shù)大于滲透實(shí)驗(yàn)結(jié)果。另外，考慮到底部礫石材料的影響，使邊坡底部滲流特性發(fā)生變化，數(shù)值分析結(jié)果與實(shí)驗(yàn)值相比大幅減小，難以取得整合性，從而對(duì)模型邊坡整體的滲透系數(shù)產(chǎn)生較大影響。圖10為坡面發(fā)生變形時(shí)，通過滲透系數(shù)值分析得到的流速矢量。

模型降雨過程中，在邊坡發(fā)生變形區(qū)域附近，流速矢量在邊坡表面變形位置的下方開始變大，且流速矢量沿水平向外變化，說明邊坡底部發(fā)生初期變形可能是受滲透壓作用的影響。

3彈塑性有限元對(duì)初期變形的再現(xiàn)模擬

3.1數(shù)值模型和條件

圖11為數(shù)值模型截面的網(wǎng)格視圖以及邊界條件。

由圖可知，邊界條件是左右兩端的垂直方向?yàn)樽杂啥?，水平方向?yàn)楣潭ǘ?。此外，分析條件將邊坡水位設(shè)定為初期變形時(shí)從壓力計(jì)上得到的數(shù)值。在此水位狀態(tài)下，先進(jìn)行自重應(yīng)力分析，然后根據(jù)強(qiáng)度折減法確定安全系數(shù)和變形行為。巖土材料的本構(gòu)方程遵從Mohr-Coulomb準(zhǔn)則。非飽和土的有效應(yīng)力根據(jù)Bishop[11]提出：

σ′=（σ-ua）+xg（ua-uw）=

σnet+xgs，（4）

式中：σ′為有效應(yīng)力;σ為全應(yīng)力;ua為孔隙氣壓;uw為孔隙水壓;s為吸力;x為Bishop有效應(yīng)力系數(shù)，0

數(shù)值模擬中使用的物理參數(shù)如表2所示。

邊坡發(fā)生初期變形前后的狀態(tài)如圖12所示。

可以觀察到邊坡底部設(shè)置的礫石向末端方向移動(dòng)。此時(shí)坡面表層的位移量非常小，水平方向的位移量為0.5～1.5 mm。在本研究中，通過考慮微應(yīng)變到小應(yīng)變區(qū)域的變形系數(shù)來再現(xiàn)圖12所示A點(diǎn)的位移實(shí)驗(yàn)結(jié)果。

巖土材料的變形系數(shù)隨著有效應(yīng)力的減小和應(yīng)變的增大而減小并表現(xiàn)出強(qiáng)烈的非線性[12-13]。變形系數(shù)的大小對(duì)應(yīng)力和應(yīng)變分析很重要，因?yàn)闀?huì)影響初期模擬計(jì)算。本研究中數(shù)值模型的微應(yīng)變至小應(yīng)變的變形系數(shù)可通過常規(guī)三軸壓縮實(shí)驗(yàn)獲得。本研究采用經(jīng)驗(yàn)公式（6）和公式（7）[14]計(jì)算了邊坡模型和支撐地層彈性區(qū)域的變形系數(shù)。

Gmax=1 400×f（e）×σ′m0.4，（6）

f（e）=（2.17-e）2/（1+e），（7）

式中：Gmax為微應(yīng)變區(qū)的剪切模量;f（e）為孔隙比函數(shù);σ′m為平均有效應(yīng)力;e為孔隙比，由邊坡和下部支撐土層的ρd和ρs計(jì)算得出。

3.2數(shù)值模擬結(jié)果分析

剪切模量變化差值ΔGmax與坡面變形量的關(guān)系如圖13所示。

圖10和圖12中坡面變形，A點(diǎn)產(chǎn)生的水平位移為ux。在0.7Gmax左右時(shí)實(shí)驗(yàn)與數(shù)值模擬結(jié)果基本一致，表示可以通過降低無水位狀態(tài)時(shí)的Gmax值來預(yù)測(cè)邊坡初期變形狀態(tài)。非飽和狀態(tài)下，Gmax的增減除了受滲透壓的影響外，也需考慮由于飽和度增加吸力降低導(dǎo)致有效應(yīng)力減少的可能性[15]。

采用0.7Gmax得到的邊坡底部變形的數(shù)值分析結(jié)果如圖14所示。

數(shù)值模擬的邊坡底部發(fā)生變形，表明數(shù)值分析結(jié)果與實(shí)驗(yàn)結(jié)果一致。在上述模擬中，假設(shè)變形系數(shù)減小的區(qū)域是整個(gè)模型邊坡，為了確定模型的平均變形系數(shù)，式（6）中的垂直有效應(yīng)力取值為模型中心所受垂直應(yīng)力（H=300 mm）。

水平位移分布的模擬結(jié)果如圖15所示。發(fā)生初期變形后持續(xù)降雨，650 min后發(fā)生大面積崩塌。發(fā)生大面積崩塌時(shí)模型內(nèi)水位相比圖15的水位仍處于上升狀態(tài)。從圖15可以看出，模擬結(jié)果發(fā)現(xiàn)發(fā)生大變形的區(qū)域與圖16所示實(shí)驗(yàn)中的最終破壞區(qū)域基本一致。因此，在初期變形發(fā)生時(shí)，邊坡變形系數(shù)（剛度）下降引起的微小變形極有可能形成潛在的滑移面，從而導(dǎo)致大面積崩塌。邊坡出現(xiàn)初始變形后進(jìn)一步降雨，地下水位上升引起的滲透壓作用使得有效應(yīng)力減小，從而導(dǎo)致邊坡逐漸失穩(wěn)。為了防止土體結(jié)構(gòu)逐漸失穩(wěn)而導(dǎo)致大面積崩塌，建議采取坡體排水等措施來防止邊坡發(fā)生初期變形。由于邊坡的初始微小變形很可能與潛在滑動(dòng)面有關(guān)，因此，如果檢測(cè)到這種微小變形，就必須考慮邊坡發(fā)生崩塌或滑坡的可能性。

4結(jié)論

基于粉質(zhì)砂土模型邊坡降雨實(shí)驗(yàn)，通過土體飽和度和孔隙水壓力的變化分析了邊坡失穩(wěn)的過程，通過飽和及非飽和滲流模擬再現(xiàn)邊坡水位變化過程，采用彈塑性有限元法再現(xiàn)邊坡初始變形狀態(tài)，具體結(jié)論如下。

1）通過對(duì)滲流過程中地下水位上升進(jìn)行分析，并使用彈塑性有限元法對(duì)模型實(shí)驗(yàn)隨水位上升的初期變形進(jìn)行數(shù)值模擬，發(fā)現(xiàn)：必須使用比滲透實(shí)驗(yàn)得到的滲透系數(shù)更大的數(shù)值進(jìn)行計(jì)算，才能在數(shù)值分析中再現(xiàn)降雨過程中的邊坡水位變化情況。

2）在模擬邊坡變形系數(shù)降低的彈塑性有限元分析中，通過將彈性區(qū)域的變形系數(shù)降低30%得到與實(shí)驗(yàn)值相一致的數(shù)值解，由此可以看出，可以采用降低變形系數(shù)的數(shù)值模擬方法來預(yù)測(cè)邊坡初期變形狀態(tài)。

3）當(dāng)前研究結(jié)果僅限于模型實(shí)驗(yàn)的再現(xiàn)且數(shù)值模擬基于邊坡整體變形系數(shù)降低，而實(shí)際邊坡破壞時(shí)，在滲透壓作用下只有邊坡底部局部變形系數(shù)的降低。

在今后研究中，需要對(duì)降雨實(shí)驗(yàn)和實(shí)際因降雨引起滑坡的案例進(jìn)行更詳細(xì)的數(shù)值分析，在此基礎(chǔ)上建立更準(zhǔn)確的邊坡崩塌預(yù)測(cè)方法，提出更有效的滑坡防治措施。

參考文獻(xiàn)/References：

[1]ORENSE R P.Slope failures triggered by heavy rainfall[J].Philippine Engineering Journal，2004，25（2）：73-90.

[2]TOHARI A，NISHIGAKI M，KOMATSU M.Laboratory rainfall-induced slope failure with moisture content measurement[J].Journal of Geotechnical and Geoenvironmental Engineering，2007，133（5）：575-587.

[3]張明禮，溫智，董建華，等.考慮降雨作用的多年凍土區(qū)不同地表土質(zhì)活動(dòng)層水熱過程差異分析[J].巖土力學(xué)，2020，41（5）：1549-1559.

ZHANG Mingli，WEN Zhi，DONG Jianhua，et al.Response of hydrothermal activity in different types of soil at ground surface to rainfall in permafrost region[J].Rock and Soil Mechanics，2020，41（5）：1549-1559.

[4]羅中寶.液壓離心地震模擬振動(dòng)臺(tái)模型的建立及其控制策略研究[D].哈爾濱：哈爾濱工業(yè)大學(xué)，2015.

LUO Zhongbao.Modelling and Control Strategr Study of Hydraulic Centrifuge Shaker Used in Earthquake Simulation[D]Harbin：Harbin Institute of Technology，2015.

[5]LING H，LING H I.Centrifuge model simulations of rainfall-induced slope instability[J].Journal of Geotechnical and Geoenvironmental Engineering，2012，138（9）：1151-1157.

[6]胡嫚，吳飛，汪時(shí)機(jī)，等.基于光滑粒子法邊坡失穩(wěn)影響范圍的模擬[J].重慶大學(xué)學(xué)報(bào)，2019，42（5）：56-65.

HU Man，WU Fei，WANG Shiji，et al.Modeling of influenced area after slope failure based on smoothed particle hydrodynamics （SPH）[J].Journal of Chongqing University，2019，42（5）：56-65.

[7]包一丁.排土場(chǎng)邊坡穩(wěn)定性評(píng)價(jià)及潛在滑坡災(zāi)害范圍預(yù)測(cè)——以朱家包包排土場(chǎng)為例[D].長(zhǎng)春：吉林大學(xué)，2020.

BAO Yiding.Stability Analysis and Landslide Failure Potential Forecasting of the Waste Dump：Zhujiabaobao Waste Dump Case Study[D]Changchun：Jilin University，2020.

[8]BUI H H，F(xiàn)UKAGAWA R，SAKO K，et al.Slope stability analysis and discontinuous slope failure simulation by elasto-plastic smoothed particle hydrodynamics（SPH）[J].Géotechnique，2011，61（7）：565-574.

[9]VAN GENUCHTEN M T.A closed form equation for predicting the hydraulic conductivity of unsaturated soils[J].Soil Science Society of America Journal，1980，44（5）：892-898.

[10]MUALEM Y.A new model for predicting the hydraulic conductivity of unsaturated porous media[J].Water Resources Research，1976，12（3）：513-522.

[11]BISHOP A W，BLIGHT G E.Some aspects of effective stress in saturated and partly saturated soils[J].Geotechnique，1963，13（3）：177-197.

[12]TATSUOKA F.Laboratory stress-strain tests for developments in geotechnical engineering research and practice[C]//Proceedings of the Fifth International Symposium on Deformation Characteristics of Geomaterials.Amsterdam：IOS Press，2011：3-50.

[13]TATSUOKA F，IWASAKI T，F(xiàn)UKUSHIMA S，et al.Stress conditions and stress histories affecting shear modulus and damping of sand under cyclic loading[J].Soils and Foundations，1979，19（2）：29-43.

[14]IWASAKI T，TATSUOKA F.Effects of grain size and grading on dynamic shear moduli of sands[J].Soils and Foundations，1977，17（3）：19-35.

[15]KAWAJIRI S，SHIBUYA S，KATO S，et al.Effects of matric suction on elastic shear modulus for three unsaturated soils[C]//Proceedings of the 5th International Conference on Unsaturated Soils.[S.l.]：[S.n.]，2010：271-275.