王　俊,林國進,唐　協(xié),何　川

(1.西南交通大學交通隧道工程教育部重點實驗室,四川 成都 610031; 2.四川省交通運輸廳公路規(guī)劃勘察設計研究院,四川 成都 610041)

砂土地層是我國城市盾構隧道建設中遇到的一種典型地質條件,具有離散性強、黏聚力弱等特點,對外界擾動反應靈敏.盾構在該類介質中掘進時若支護壓力不足可能誘發(fā)掌子面坍塌,給國民經濟造成重大損失.因此開展砂土地層盾構掘進掌子面的穩(wěn)定性研究是十分必要的.

目前,采用理論分析、模型試驗以及數值模擬等方法對盾構隧道掌子面穩(wěn)定性進行了研究.理論分析包括極限平衡理論與極限分析理論,前者需預先假定失穩(wěn)區(qū)大小與形狀,如對數螺旋曲線[1]、半球體和半圓[2]、楔形體[3-5]等.后者則需要假設應力場分布[6-8](下限解)或者速度場分布[8-9](上限解).

模型試驗主要通過應力[10]或應變[11-13]的控制方式研究隧道掌子面穩(wěn)定性.Chambon和Corté[10]通過逐步減小支護壓力并結合地層變形情況得到掌子面極限支護壓力;Kirsch與Idinger等[11-12]則通過后撤擋板模擬盾構開挖,同時監(jiān)控作用在擋板上的土壓力,并根據開挖狀態(tài)將其視為極限或者殘余支護壓力.

隨著計算機技術的不斷發(fā)展,數值模擬在隧道掌子面穩(wěn)定性研究中得到了廣泛運用.有限元方面,Vermeer等[14]認為排水條件下極限支護壓力隨內摩擦角的增大而減小;Li等[15]采用FLAC3D研究了軟土地層大直徑泥水盾構掘進掌子面穩(wěn)定性問題.離散元方面[16-20],王俊等[17]采用PFC2D研究了砂卵石地層中泥水盾構掌子面穩(wěn)定性;Chen等[19]通過后側擋板模擬了砂土地層盾構隧道掌子面失穩(wěn).

上述研究雖取得諸多成果,但仍有一定不足:如理論分析主要適用于簡單本構土體,一般不考慮土拱效應與隧道開挖過程的影響;模型試驗能從宏觀上把握掌子面失穩(wěn)現象,但限于觀測手段難以解釋失穩(wěn)機理且不易開展參數研究;數值模擬主要采用有限元方法開展,計算結果受本構模型選取、網格劃分精度等因素影響較大,且在模擬土體大變形問題上存在困難.既有的三維離散元研究也存在局限性:如Chen等[19]采用的研究方法難以模擬隧道開挖過程,且不能直接得到掌子面極限支護壓力;Melis Maynar和Medina Rodrigue[20]則采用地層黏聚力評價掌子面穩(wěn)定性,研究成果不能用于砂土地層.

鑒于此,本研究以Chambon和Corté[10]開展的室內模型試驗為基礎,采用三維離散元方法建立了與之相匹配的數值模型,從細觀角度解釋了砂土中盾構隧道開挖面失穩(wěn)機理,研究了極限支護壓力、失穩(wěn)區(qū)分布、土拱效應和地層移動規(guī)律.數值模型中引入了三維動態(tài)柔性應力邊界,該應力邊界可有效模擬室內試驗中倉內介質(空氣、流體)對開挖面的支撐效應,將其抽象為作用在掌子面顆粒上的指定支護壓力,模擬過程中分步等量減小該壓力直到隧道失穩(wěn).與既有研究相比,本文通過刪除進入隧道內部的土顆粒模擬隧道開挖,考慮了該施工行為對掌子面穩(wěn)定性的影響,采用的三維柔性應力邊界可精確得到極限支護壓力.

1　三維動態(tài)柔性應力邊界

1.1　室內試驗簡介

Chambon與Corté[10]開展的模型試驗如圖1所示.模型箱尺寸為1.20 m(長)×0.80 m(寬)×0.36 m(高),采用Fontainebleau砂進行試驗,主要物理力學性質為:顆粒平均粒徑d50=0.17 mm;不均勻系數Cu=1.47;重度為15.3～16.1 kN/m3;內摩擦角為38°～42°;黏聚力為0～5 kPa.采用靜水壓力或空氣壓力模擬試驗中掌子面處的支護壓力,考慮了C/D=0.5,1.0,2.0,4.0 等4種埋深工況,其中:C為隧道埋深,D=0.1 m 為隧道直徑.試驗考慮了50g、100g和130g3種加速度情況,對應的原型隧道直徑分別為5、10 m 和 13 m.

室內試驗中首先將掌子面處的支護壓力設為隧道軸線處的主動土壓力,此后逐步減少該壓力,模擬隧道失穩(wěn)結合地層變形情況得到極限支護壓力與失穩(wěn)區(qū)分布.

1.2　三維動態(tài)柔性邊界算法

采用三維離散元方法PFC3D建立Chambon和Corté室內試驗數值模型的主要難點是如何描述模型試驗中靜水壓力或空氣壓力對掌子面的支護效應.由于水和空氣均是流動性較強的介質,不能承受剪切力,故其對隧道的支護效應可視為作用在掌子面上的法向力,倉內介質與開挖面的相互作用可簡化為三維應力邊界問題,且該應力邊界隨隧道掌子面處土體位移變形而不斷演化.

圖1　Chambon和Corté室內模型試驗Fig.1　Sketch of Chambon and Corté’s model test

為在PFC3D模型中構建該應力邊界,通過內置Fish語言[21]進行編程將倉內介質對掌子面的支護效應抽象成作用在掌子面顆粒上的指定法向壓力.如圖2(a)所示,程序實現的基本思想是假設在盾構正后方存在一個邊長與盾構直徑相等的平面,并對其進行離散化處理,得到一系列晶格,所有的顆粒均沿隧道軸線向離散晶格投影,得到該方向上距對應晶格最近的所有顆粒(即掌子面顆粒),再根據顆粒在對應晶格上的投影面積Ap(圖2(b))給掌子面顆粒施加指定的支護力f=pAp(p為倉內支護壓力,kPa),便實現了給掌子面施加支護,構建了倉內介質與掌子面土體之間相互作用的應力邊界.為了便于理解,可以將晶格視作向掌子面發(fā)射平行光的光源,所有能接收光照的顆粒就是掌子面顆粒,其所受的支護力與光照面積成正比.

劃分晶格時需注意砂土顆粒平均直徑(davg)與晶格尺寸(l)之間的大小關系,davg/l過小則計算誤差大,過大則計算耗時,經過多次試算最終選取davg/l=16.掌子面顆粒與離散晶格的相對關系如圖3所示.

(a)離散網格示意圖(b)土壓力施加圖2　土倉壓力建立算法示意圖Fig．2　Sketchoftheearthpressureinstallation

圖3　掌子面顆粒與離散晶格相對關系Fig.3　Relationship between grid cell and sand particle

2　三維離散元數值模型

2.1　細觀參數標定

采用三軸試驗對模型中采用的砂土顆粒進行了標定.考慮到隧道埋深,標定圍壓σ3選為50、100、200 kPa.最終選用的細觀力學參數如表1所示,標定后的宏觀力學參數內摩擦角φ=38°,黏聚力c=0 kPa,與室內模型試驗采用的Fontainebleau砂的物理力學性質基本一致.1#材料標定過程的應力應變曲線如圖4所示.

表1　PFC3D細觀力學參數Tab.1　Calibrated PFC3D microscopic parameters

圖4　材料細觀參數標定Fig.4　Calibration of microscopic parameters

本次三軸試驗對傳統(tǒng)標定方法進行了改進,采用前述三維柔性應力邊界代替剛性wall單元構建了側向約束,有效模擬了室內試驗中橡皮膜對試樣的柔性包裹.如圖5(a)所示,假設在試樣周圍存在兩個正交平面,并對其進行離散化處理,得到一系列晶格.所有土顆粒均沿x和z方向向晶格i(1≤i≤n,n為晶格總數)投影,得到距離晶格i最近和最遠的顆粒.當確定了距任意晶格最遠和最近的土顆粒便得到了試樣中的邊界顆粒.作用在邊界顆粒上的力則為指定圍壓與顆粒橫截面積的乘積,再根據顆粒形心坐標與試樣中心的夾角將該力沿x與z方向進行分解并施加在對應顆粒上,便實現給所有邊界顆粒施加相等的徑向壓力.隨著試樣不斷加載,邊界顆粒的位置會發(fā)生變化,數值三軸試驗中每隔100步執(zhí)行一次三維柔性應力邊界算法,更新邊界顆粒并調整作用在其上的圍壓.

最終剪切破壞的試樣如圖5(b)所示.

為提高計算效率,在隧道軸線范圍內的地層由1#土顆粒構成,其余范圍內的地層由2#土顆粒組成,兩種顆粒均采用線性接觸模型,顆粒粒徑服從 0～1均勻分布.

(a)網格劃分示意圖(b)剪切破壞試樣圖5　三軸試驗示意圖Fig．5　Sketchoftriaxialtest

2.2　PFC3D數值模型

如圖6所示,綜合考慮邊界效應與計算效率,數值模型尺寸為:長(L)×寬(B)=18 m×25 m,高度(H)由隧道埋深確定,共考慮了C/D=0.5、1.0、2.0、4.0 4種埋深工況,隧道直徑D=5.0 m.隧道軸線上下方各3.5 m范圍內(粒徑滿足8～12cm的顆粒分布范圍m=7.0 m)的砂土顆粒粒徑分布為8～12 cm,其余部位為15～20 cm.盾構埋入地層中的長度S1為7.0 m,地層外的長度S2為2.5 m,隧道下方土體厚度W為4 m.為了得到隧道失穩(wěn)過程中地層變形情況,模型中布設了如圖6所示的位移測點.

數值模擬按以下步驟進行:

步驟1采用wall單元生成地層模型邊界后(wall單元法向與切向剛度均為1×108N/m,摩擦因數為0),根據標定的細觀參數,采用“落雨法”分層生成地層模型,在重力作用下完成固結后刪除地層頂部少量浮皮顆粒;

步驟2在指定范圍內(y≤0,x2+z2≤2.52)刪除地層顆粒并采用wall單元生成盾構機,其法向與切向剛度均為1×1011N/m,摩擦因數為0.2;

步驟3利用三維動態(tài)柔性邊界算法識別掌子面顆粒,并施加與隧道軸線處靜止土壓力相等的初始支護壓力,此后分步等量減少該支護壓力直至0 kPa,同時密切監(jiān)控掌子面前方土體位移變化,當測點位移突然增大時,可以認為此時的支護壓力為掌子面極限支護壓力[10].

(a)縱剖面(b)縱剖面圖6　模型尺寸與測點布置示意圖Fig．6　Sketchofmodelandmonitoringpointdistribution

為提高計算精確度,每一級支護壓力下均迭代50 000步,保證模型達到平衡狀態(tài).C/D=2.0時的數值模型見圖7.

圖7　PFC3D數值模型Fig.7　PFC3D numerical model

3　數值模擬結果分析與討論

通過分析數值計算結果中地層變形以及土體應力場數據,得到了掌子面極限支護力、掌子面失穩(wěn)破壞模式、地層移動與掌子面附近土體應力變化規(guī)律.本文主要介紹開挖面極限支護力、失穩(wěn)區(qū)范圍與開挖面前方土拱效應發(fā)展規(guī)律等.

3.1　數值模型合理性驗證

將數值計算結果與室內模型試驗相對比驗證其合理性.為了得到極限支護壓力,在減小支護壓力的同時密切記錄地層變形情況,得到兩者之間的關系曲線.限于篇幅本文只給出了C/D=2.0時地中測點位移隨支護壓力減小變化情況(圖8).該測點位于隧道拱頂正前方1.0 m,坐標為(0.0,1.0,2.5),由文獻[10-12]知,該點位于失穩(wěn)區(qū)中心位置,能及時反應隧道失穩(wěn).

圖8　地中測點沉降與支護壓力關系曲線Fig.8　Relationship between subsurface settlement and support pressure

從圖8可以看出,測點位移隨支護壓力p變化曲線可以分為兩個階段:初始階段測點位移隨p減小逐漸增大,但變化幅度很小;當p減小到特定值時,測點位移突然快速增大,由文獻[10]可判定p為極限支護壓力(pf).按此辦法,本研究給出了所有工況的掌子面極限支護壓力,具體結果如表2所示.

表2　不同埋深條件下掌子面極限支護壓力Tab.2　Limit support pressure under various C/D

如圖9所示,引入無量綱系數pf/(γD)將離散元計算與Chambon和Corté[10]模型試驗結果進行了比較,γ為土體重度.

由圖9可知,兩種研究中極限支護壓力隨隧道埋深變化規(guī)律是相同的,可以分為兩個階段:C/D≤1.0時極限支護壓力隨C/D增長較快,而當C/D≥2.0時極限支護壓力趨于穩(wěn)定,基本不隨埋深變化而變化.但數值計算得到的結果絕對值大于模型試驗,出現該結果的可能原因是Chambon和Corté[10]離心試驗采用的砂土物理力學性質具有一定的不確定性,內摩擦角在38°～42°范圍內變化大于數值模型中砂土的內摩擦角(φ=38°),還可能具有一定的黏聚力(0≤c≤5 kPa),而數值模型中砂土黏聚力為0.由Vermeer等[14]的研究可知掌子面極限支護壓力隨砂土地層的內摩擦角減小而增大,此外由文獻[7-9]知,黏聚力可以有效提高地層自穩(wěn)性,減小極限支護壓力,故數值計算與模型試驗之間的差異是合理的.

圖9　數值計算結果與其他研究比較Fig.9　Comparison between numerical results and other researchs

此外,將離散元計算結果與理論研究[4,8-9]和Kirsch模型試驗[11]結果也進行了對比.離散元計算結果略大于Leca和Dormieux[8]與Mollon等[9]給出的上限解,由極限分析中上限解的定義知該結論是合理的;離散元計算結果與Anagnostou和Kovari[4]提出的三維楔形體模型吻合良好,當C/D≤1.0時兩者結果基本相同,當C/D=2.0時前者略小于后者,出現該現象的可能原因是當C/D≥2.0時隧道拱頂上方出現土拱效應,提高了隧道的自穩(wěn)能力,而楔形體模型中由于沒有考慮該效應,掌子面極限支護壓力大于數值計算結果,關于土拱效應的具體分析可見3.5小節(jié).Kirsch[11]模型試驗中采用的砂土摩擦角為32.5°,小于本模型中的38.0°,Vermeer等[14]認為掌子面極限支護壓力隨砂土地層的內摩擦角減小而增大,故離散元結果小于該模型試驗也是合理的.綜上所述,研究采用的數值模型是合理可行的.

定義支護壓力比η=p/p0,其中,p0為初始狀態(tài)隧道軸線處靜止土壓力,則4種工況下極限支護壓力比ηf分別為12.0%、11.7%、9.6%、5.8%,ηf隨C/D的增加而減小.

3.2　失穩(wěn)區(qū)分布

明確失穩(wěn)區(qū)形狀與大小分布是掌子面穩(wěn)定性研究中的重點,特別是在極限平衡理論分析中.本文通過以下兩個步驟得到失穩(wěn)區(qū)分布:

步驟1記錄圖6中所有測點的豎向位移;

步驟2根據同一斷面處測點位移梯度判定失穩(wěn)區(qū)范圍,依次連線所有位移突增點得到失穩(wěn)區(qū)分布.

圖10給出了隧道失穩(wěn)時典型斷面處(C/D=2.0,y=0.5,z=3.0)測點位移分布情況.

圖10　典型橫斷面上地層變形情況Fig.10　Typical vertical displacement curve at transverse profile

如圖11所示,采用上述辦法給出了4種埋深工況下隧道失穩(wěn)區(qū)分布形態(tài).

由圖11可知,對C/D=0.5來講,p=pf,0 kPa時失穩(wěn)區(qū)均呈煙囪狀,且發(fā)展到地表,數值計算得到的分布范圍略大于模型試驗,出現該現象的可能原因是模型試驗中砂土內摩擦角大于數值模型,且可能存在一定的黏聚力.

計算結果表明:失穩(wěn)區(qū)從拱底開始向上發(fā)展,p=pf時發(fā)展至掌子面前方約0.77D處,p=0 kPa時發(fā)展至1.02D,此外失穩(wěn)區(qū)從拱頂處還向隧道后方延伸一定范圍;C/D=1.0時,失穩(wěn)區(qū)同樣發(fā)展至地表,與模型試驗結果不同,具體形態(tài)與C/D=0.5基本相同,但分布范圍更大;C/D=2.0,4.0時,失穩(wěn)區(qū)沒有發(fā)展到地表,分布范圍略大于模型試驗,兩種埋深條件下失穩(wěn)區(qū)均在拱頂后方有一定范圍的延伸.

(a)C/D=0．5(b)C/D=1．0(c)C/D=2．0(d)C/D=4．0圖11　失穩(wěn)區(qū)縱剖面圖Fig．11　Longitudinalprofileoffailurezone

3.3　地表沉降

圖12給出了C/D=2.0時地表典型測點(坐標為(0.0,1.0,12.5))沉降隨支護壓力變化曲線.

由圖12可知,相比于地中測點,地表測點沉降具有滯后性,地表沉降突增點對應的支護壓力pk

圖12　地表測點沉降與支護壓力關系曲線Fig.12　Relationship between surface settlement and support pressure

為了整體把握地表沉降情況,給出了4種工況下沉降三維曲面.如圖13所示,地表沉降范圍隨埋深增大而增大,最大沉降均出現在開挖面前方約0.16D(y=0.8 m)處,且其絕對值隨埋深增大減小.

表3　不同埋深條件下地表位移突增時的支護壓力Tab.3　pk under various C/D

由圖13(a)可知:由于隧道埋深較淺,當C/D=0.5時,失穩(wěn)區(qū)基本呈煙囪狀發(fā)展至地表,地表沉降主要集中在隧道軸線兩側各4 m (0.8D)以及前方5 m (1.0D)范圍內,同時向開挖面后方發(fā)展了大約2 m (0.4D),最大沉降達到了5.85 m;C/D=1.0時,地表沉降分布范圍集中在軸線兩側各5 m (1.0D)以及前方6 m (1.2D)范圍內,同時向隧道開挖面后方發(fā)展了大約3 m (0.6D),地表最大沉降為2.42 m;C/D=2.0時,地表沉降形態(tài)與C/D=1.0比較相似,但沉降區(qū)范圍更廣,地表最大沉降為0.41 m;C/D=4.0時,地表沒有出現明顯的規(guī)律.

(a)C/D=0．5(b)C/D=1．0(c)C/D=2．0(d)C/D=4．0圖13　地表沉降等高線圖Fig．13　Surfacesettlementcontour

3.4　地中沉降

采用地中沉降槽寬度參數K研究了不同埋深條件下地中沉降情況,K=i/(z0-z),其中:i為沉降槽寬度系數;z和z0分別為地中某一水平面與隧道中心距離地表的高度.

圖14為地中沉降槽寬度參數隨深度比的變化情況,由圖14可知:砂土地層的地中沉降槽寬度參數隨深度比的變化規(guī)律與既有針對黏土和復合地層的研究結果差別較大,隨著深度比的增大,C/D=0.5,1.0時地中沉降槽寬度參數基本呈線性增加,相同深度比條件下,C/D=0.5的地中沉降槽寬度參數要小于C/D=1.0,且兩種地層的地中沉降槽寬度參數均大于復合地層地中沉降槽寬度參數,同時又都小于黏土地層地中沉降槽寬度參數;C/D=2.0時,可能由于拱效應,沉降槽寬度參數變化沒有明顯規(guī)律,但仍基本處于復合地層與黏土地層的包絡線內.

圖14　地中沉降槽寬度參數隨深度比的變化Fig.14　Width parameter of subsurface settlement trough versus buried depth ratio

3.5　土拱效應分析

當C/D=2.0，4.0時,在掌子面上前方出現了土拱效應,有效地提高了隧道穩(wěn)定性.為了定量解釋該現象,圖15、16中給出了兩種工況下隧道掌子面前方0.5D處p=p0，pf，0 kPa時,垂直隧道軸線方向水平應力σxx、沿隧道軸向方向水平應力σyy、豎向應力σzz沿隧道埋深分布情況.

由圖15可知,C/D=2.0時,應力分布沿埋深方向大致分為3個階段.以σyy為例:(1) 由于隧道開挖,從隧道拱底(A點)至隧道拱頂(B點)應力釋放明顯,輪廓范圍內應力不同程度減小;(2) 從拱頂至隧道上方一定范圍內應力逐漸增加,并在拱頂上方約 4 m處達到極值;(3) 從C點至地表應力逐漸降為0.

(a)σxx(b)σyy(c)σzz圖15　C/D=2．0掌子面前方應力沿隧道埋深重分布(y=0．5D)Fig．15　StressredistributionatC/D=2．0alongburieddepthinfrontoftunnelface(y=0．5D)

(a)σxx(b)σyy(c)σzz圖16　C/D=4．0掌子面前方應力沿隧道埋深重分布(y=0．5D)Fig．16　StressredistributionatC/D=4．0alongburieddepthinfrontoftunnelface(y=0．5D)

由Chen等[13]研究可知,可通過定義水平應力集中系數λ判定土拱分布范圍,λ=σ/σ0,其中:σ為任意時刻土壓力;σ0為初始土壓力.由于卸載效應,拱頂上方失穩(wěn)區(qū)內有λ<1.0,而土拱將承擔失穩(wěn)區(qū)與兩側土體轉移來的荷載,故有λ>1.0.

由圖15可知:p=pf時土拱存在于拱頂上方約0.7D～1.3D范圍內;當p從pf降至0 kPa時,土拱有向地表發(fā)展的趨勢,與Chen等[13]室內試驗結果基本吻合.

由圖16可知:對C/D=4.0工況來講,p=pf時主應力沿深度變化規(guī)律與C/D=2.0工況基本一致,塌落拱出現在隧道上方約0.9D～2.3D范圍內;p降至 0 kPa后,土拱仍然向地表發(fā)展了一段距離.

4　結　論

以Chambon和Corté室內試驗為背景,采用PFC3D開展離散元數值模擬從細觀角度解釋了砂土地層盾構隧道失穩(wěn)機理.通過與理論分析和室內試驗進行對比驗證了數值模型的合理性.主要結論如下:

(1) 掌子面極限支護壓力pf與隧道埋深相關,C/D≤1.0時,極限支護壓力隨C/D增長較快;C/D≥2.0時,極限支護壓力趨于穩(wěn)定,基本不隨埋深變化.極限支護壓力比隨埋深增加而減少.

(2) 數值計算得到的失穩(wěn)區(qū)分布范圍稍大于模型試驗,且失穩(wěn)區(qū)均向隧道后方發(fā)展一段距離.C/D=0.5,1.0時,掌子面失穩(wěn)向上發(fā)展引起地表坍塌;C/D=2.0,4.0時,在拱頂上方形成了穩(wěn)定的塌落拱,塌落拱高度分別約為0.7D～1.3D與0.9D～2.3D.

(3) 地表沉降突增點對應的支護壓力pk

(4) 隧道埋深越淺,地表沉降范圍越小,最大沉降值越大,最大地表沉降出現在掌子面前方0.16D處.C/D=0.5,1.0時砂土地層的地中沉降槽寬度參數大于復合地層但小于黏土地層,且隨埋深比增加線性增長.

致謝:2014年度西南交通大學博士研究生創(chuàng)新基金;第四屆西南交通大學軌道交通行業(yè)拔尖創(chuàng)新人才培育項目.

參考文獻:

[1]MURAYAMA S,ENDO M,HASHIBA T,et al.Geotechnical aspects for the excavating performance of the shield machines[C]∥The 21st Annual Lecture in Meeting of Japan Society of Civil Engineers.Tokyo:[s.n.],1966:134-140.

[2]KRAUSE T.Schildvortrieb mit flüsigkeits-und erdgestüzter ortsbrust[D].Brunswick:Technical University Carolo Wilhelmina,1987.

[3]HORN N.Horizontal earth pressure on the vertical surfaces of the tunnel tubes[C]∥National Conference of the Hungarian Civil Engineering Industry.Budapest:[s.n],1961:7-16.

[5]BROERE W.Tunnel face stability and new CPT application[D].Delft:Delft University,2001.

[6]ATKINSON J H,POTTS D M.Stability of a shallow circular tunnel in cohesionless soil[J].Geotechnique,1977,27(2):203-215.

[7]DAVIS E H,GUNN M J,MAIR F R,et al.The stability of shallow tunnels and underground openings in cohesive material[J].Geotechnique,1980,30(4):397-416.

[8]LECA E,DORMIEUX L.Upper and lower bound solutions for the face stability of shallow circular tunnels in frictional material[J].Geotechnique,1990,40(4):581-606.

[9]MOLLON G,DIAS D,SOUBRA A H.Face stability analysis of circular tunnels driven by a pressurized shield[J].Journal of Geotechnical and Geoenvironmental Engineering,2010,136(1):215-229.

[10]CHAMBON P,CORTé J F.Sallow tunnels in cohesionless soil:stability of tunnel face[J].Journal of Geotechnical Engineering,1994,120(7):1148-1164.

[11]KIRSCH A.Experimental investigation of the face stability of shallow tunnels in sand[J].Acta Geotechnica,2010,5(1):43-62.

[12]IDINGER G,AKLIK P,WU W,et al.Centrifuge model test on the face stability of shallow tunnel[J].Acta Geotechnica,2011,6(2):43-62.

[13]CHEN R P,LI J,KONG L G,et al.Experimental study on face stability of shield tunnel in sand[J].Tunnelling and Underground Space Technology,2013,33(1):12-21.

[14]VERMEER P A,RUSE N M,MARCHER T.Tunnel heading stability in drained ground[J].Felsbau,2002,20(6):8-18.

[15]LI Y,EMERIAULTB F,KASTNERB R,et al.Stability analysis of large slurry shield-driven tunnel in soft clay[J].Tunnelling and Underground Space Technology,2009,24(4):472-481.

[16]CUNDALL P A,STRACK O D.A discrete numerical model for granular assemblies[J].Géechnique,1979,29(1):47-65.

[17]王俊,何川,封坤,等.砂卵石地層中大斷面泥水盾構泥膜形態(tài)研究[J].現代隧道技術,2014,51(6):108-113.

WANG Jun,HE Chuan,FENG Kun,et al.Research on the dynamic behavior of the slurry membrane of a large-section slurry shield in a sandy cobble stratum[J].Modern Tunnelling Technology,2014,51(6):108-113.

[18]繆林昌,王正興,石文博.砂土盾構隧道掘進開挖面穩(wěn)定理論與顆粒流模擬研究[J].巖土工程學報,2015,37(1):98-104.

MIU Lingchang,WANG Zhengxing,SHI Wenbo.Theoretical and numerical simulations of face stability around shield tunnels in sand[J].Chinese Journal of Geotechnical Engineering,2015,37(1):98-104.

[19]CHEN R P,TANG L J,LING D S,et al.Face stability analysis of shallow shield tunnels in dry sandy ground using the discrete element method[J].Computers and Geotechnics,2011,38(2):187-195.

[20]MELIS MAYNAR M J,MEDINA RODRIGUEZ L E.Discrete numerical model for analysis of earth pressure balance tunnel excavation[J].Journal of Geotechnical and Geoenvironmental Engineering,2005,131(10):1234-1242.

[21]Itasca Consulting Group Inc.PFC3D (particle flow code in three dimensions) version 4.0 manual[M].Minneapolis:Itasca Consulting Group Inc,2001:101-235.

[22]MAIR R J,TAYLOR R N,BRACEGIRDLE A.Subsurface settlement profiles above tunnels in clays[J].Gétechnique,1993,43(2):315-320.

[23]MAIR R J.Centrifugal modeling of tunnel construction insoft clay[D].Cambridge:University of Cambridge,1979.

[24]DYER M R,HUTCHINSON M T,EVANS N.Sudden valley sewer:a case history[C]∥International Symposium on Geotechnical Aspects of Underground Construction in Soft Ground.London:[s.n.],1996:671-676.