表1 离散元材料参数
Tab.1 Material parameters of DEM sample
参数数值颗粒粒径d/mm3～6颗粒密度ρ/(kg·m-3)2600颗粒间摩擦系数μ10.5颗粒与墙体间摩擦系数μ20颗粒与墙体法向刚度kn/(N·m-1)1.5e8颗粒与墙体切向刚度ks/(N·m-1)1.0e8
摘 要：为阐明土压平衡盾构挤土与排土掘进状态对周围土体的扰动机理，采用离散元方法建立土仓不排土与排土两组盾构掘进二维数值模型，比较分析了两种工况下盾构掘进时掌子面及周围地层的细观力学变化特征.研究表明：不排土工况下，盾构对周围土层有明显挤密，土体竖向位移呈椭圆形、水平位移呈圆形向前方扩散变小；排土工况下，刀盘周围土体无明显挤密，水平位移很小，竖向位移呈不规则漏斗状向地表延伸变大.不排土对周围土体应力影响大，掌子面水平应力分布均匀；排土工况下盾构前方土体应力与初始地层接近，土仓的排土，使掌子面水平应力的分布变化较大.
关键词：无水砂层；土压平衡盾构；排土；地层响应；离散元
土压平衡盾构掘进必定会对周围地层造成扰动，产生地层变形.目前地层变形研究大多假设土仓进排土为均衡的理想状态，而事实上，土仓内进排土大多时候都无法控制平衡.研究盾构土仓进排土状态对周围地层的扰动机理，既可提升对盾构隧道深层土体扰动机理的认识，也能提高对深层土体位移场、应力场的预测和控制能力.盾构掘进引起地层响应的主要研究方法有解析法[1-2]、室内模型试验法[3-4]、有限元或有限差分法[5-6]、离散元方法[7]等.解析法公式推导严格，但对假设条件过于苛刻.模型试验得出的结论可靠，但成本高，制样复杂.有限单元法可以盾构隧道施工过程各个施工参数对地层变形的影响，是研究盾构施工引起地层响应的常用方法之一，然而对于砂卵石地层，土体具有明显的离散特性，有限单元法将其当作一个连续的弹塑性体，其本构模型无法准确表达砂卵石土体的离散特性，更难以描述盾构掘进的动态大变形过程. 离散单元法克服了有限单元法的连续性假定，深入土体的细观级别.目前已有一些离散元运用到盾构隧道领域的研究实例，Maynar等[7]结合马德里地铁工程，首次采用离散元方法对盾构开挖进行模拟，探讨了盾构推进过程中的推力、扭矩与掌子面稳定问题.Chen等[8]利用PFC3D研究了盾构隧道掌子面支护压力随掌子面移动变化关系，研究了掌子面前方土体的破坏区、土体应力、孔隙率随掌子面移动变化规律.Jiang等[9]采用PFC2D程序研究了渣土改良对于土压平衡盾构的掘进过程中隧道纵向地层变形的影响，研究了渣土在不同改良程度下盾构掘进时周围土体的变形规律.Jiang等[10-11]考虑盾尾空隙条件，分析了盾构管片周围土体的应力重分布问题，并研究了不同埋深和土体损失在盾构隧道掘进过程中管片周围土体的压力拱效应.王俊等[12]采用离散元数值模拟与室内模型试验相结合的方法，研究了砂卵石地层条件下泥水盾构的泥膜形态规律.江英超等[13]采用离散元法对砂卵石地层中盾构隧道滞后沉降的发展形成过程进行模拟，分析了施工过程中各类因素对地层滞后沉降的影响，揭示了盾构隧道滞后沉降的原因.
笔者基于离散元法，从细观角度重点研究了土压平衡盾构施工过程中不排土与排土两种土仓状态对地层的响应，对比分析了土体深层位移、土体颗粒力链分布状态、掌子面的应力特征.
1.1 初始地层的生成
基于离散元商业软件PFC2D，首先生成左、右和下3面墙体边界，上侧为自由边界，初始地层采用“落雨法”生成，即首先在墙体内生成颗粒，颗粒在重力作用下堆积形成原始地层.由于土体颗粒数目和计算机能力的限制，将实际模型按20∶1的比例缩小，生成土层1.7 m×1.7 m，颗粒数目为139 414个，模型参数与文献[10]一致，见表1.
表1 离散元材料参数
Tab.1 Material parameters of DEM sample
参数数值颗粒粒径d/mm3～6颗粒密度ρ/(kg·m-3)2600颗粒间摩擦系数μ10.5颗粒与墙体间摩擦系数μ20颗粒与墙体法向刚度kn/(N·m-1)1.5e8颗粒与墙体切向刚度ks/(N·m-1)1.0e8
原始地层在重力作用下运动直到稳定，模拟自然沉降过程.地层稳定性的判断依据主要为颗粒最终的平均速度、最大速度及特征单元的位移，当平均速度和最大速度值趋近于0(v<1×10-5 m/s)，特征单元的位移达到某一值而基本不再变化时，可认定体系达到稳定[14].地层稳定后对顶面进行削土处理，初始地层应力平衡，位移清零，留下1.7 m×1.2 m土体作为初始地层，最终孔隙率为0.204.
1.2 盾构模型
采用PFC里墙体(wall)来模拟盾壳，如图1所示，盾构直径D为0.3 m，埋深为0.6 m(2倍洞径)，因研究重点在盾构前方和上方，盾构下方保留1倍洞径，前方保留4倍洞径，盾构土仓长0.09 m，这里不考虑刀盘面板，假定刀盘完全敞开(即开口率为100%).盾构的掘进通过赋予墙体一定的速度来实现，在排土区域一定的情况下，盾构掘进速度选定0.01 m/s.在土仓下部靠近隔板的螺旋输送机位置，设置一定的删土速度以模拟螺旋输送机的排土过程，使用Fish语言编制程序，确保删土区域与盾构相对位置保持不变，同时通过控制删土速率来模拟不同的螺旋输送机排土速度，此处平均删土速度约为0.038 m/s.笔者模拟了螺机不排土与排土两种工况，分析内容包括：
1)地层深层位移变化规律；
2)力链变化；
3)掌子面应力特征.
图1 离散元盾构模型
Fig.1 DEM model of shield tunnel
2.1 土体深层位移
为研究盾构在两种工况下掘进对周围土体位移的扰动机理，编制程序提取了盾构进尺0.15 m(D/2)时周围土体的位移，并绘制云图如图2～3所示.
图2 盾构进尺D/2时土体竖直位移云图
Fig.2 Contours of vertical ground deformation while shield advances up D/2
图3 盾构进尺D/2时土体水平位移云图
Fig.3 Contours of horizontal ground deformation while shield advances up D/2
如图2(a)，在不排土工况下，掌子面前方地表发生不同程度的隆起，最大隆起区域发生在盾顶到掌子面前方D/2范围内.在盾底及以下位置，土体向下移动，这主要是由于盾构推进挤压周围土体，土颗粒被向下挤压所致.土体竖直位移的扩散不管是向上隆起还是向下沉降均呈椭圆状，且都是刀盘附近量值大，距刀盘较远处量值小.如图2(b)所示，排土工况下，地表土体均发生沉降，竖向沉降值等值线呈倾斜不对称漏斗状.在盾构机身上方，竖向沉降值较小，主要由土颗粒自重引起.在掌子面前方，沉降值较大，主要因土体向土仓内移动，引起盾构前方土体损失所致.由图3(a)可以看出，在不排土工况下，刀盘前方明显存在挤土效应，且挤土范围表现为圆形扩散形态，掌子面附近水平位移大，远处水平位移小.原因是在不排土工况下，掌子面附近土体所受挤压应力大，远处所受挤压应力渐渐减小.如图3(b)所示，在排土工况下，除刀盘附近土体有一定水平扰动外，在其他区域土体的水平位移都很小，这说明土仓排土可有效地控制盾构刀盘前方土体的水平位移.
为分析盾构在掘进过程中两种工况下土体的变形情况，选定盾构开始掘进前掌子面前方1D位置作为监测断面，编制程序提取盾构推进5 cm(D/6)、10 cm(D/3)、15 cm(D/2)时此断面的水平和竖向深层位移值并作曲线，如图4～5所示.
由图4(a)可以看出，不排土工况下，掌子面前方土体隆起最大值发生在盾顶位置(埋深0.6 m处)，从盾顶到地表和盾底(埋深0.9 m)竖向位移逐渐减小，地表仍向上隆起，盾底土体变为下沉状态，盾底为沉降值最大位置.在盾构下方，沉降值逐渐变小，在盾底下方1D处，土体几乎无竖向位移.由图4(b)所示，排土工况下，由于土仓排土量大于进土量，引起盾构前方及上方土体损失，地层会发生沉降.埋深在0～0.6 m的土体，沉降量值随着距地表深度减小逐渐加大，即地表沉降量大，土体深部沉降量小，原因是沉降时上部土体的沉降存在一个叠加效应，地表沉降为各土层的累积沉降.而埋深在0.9 m以下的土体竖向位移很小，这是因为在排土工况下，土体未发生明显挤密.
图4 盾构始发掌子面前方1D截面竖直位移曲线
Fig.4 Curves of vertical displacements apart away 1D from cutter face launching
由图5(a)可以看出，不排土工况下，监测断面在盾顶之上(埋深0～0.6 m)，土体发生的水平位移很小；刀盘前方(埋深0.6～0.9 m)土体的水平位移较大；盾底下方(埋深0.9 m以下)，土体水平位移也很小.分析原因为：对于开口率较大的盾构，在挤土的时候，由于土仓的存在，会对盾构上部及下部土体水平向的稳定有一个保护作用，此时土体的水平位移主要发生在盾构刀盘前方.在盾构掘进过程中，土体的水平位移也逐渐增大，挤土效应更加明显.由图5(b)可以看出，在排土工况下，土体水平方向的位移很小，在刀盘附近，存在一定的扰动，但扰动量值小.在盾构掘进过程中，土体的水平位移几乎不变，这说明盾构土仓排土对刀盘前方土体水平位移控制效果明显.
2.2 力链对比
力链是颗粒物质力学概念，指外荷载通过颗粒接触传递力的大小及方向的路径.其中力链的长短代表颗粒间接触力的大小，方向代表接触力的传递方向，在同一力链图中，力链的长短是相对的.选取初始地层及盾构推进D/2时土体颗粒的力链分布，如图6所示，在不排土与排土情况下，力链的差别很大.不排土工况下，颗粒间接触力在刀盘附近较大，且由刀盘附近向前方、上下方延伸，并逐渐减弱；排土工况下，刀盘附近由于土体进入土仓，会发生接触力的释放，但这个接触力减弱的范围十分有限，在刀盘前方大多数位置，力链的分布与初始地层差别不大，这说明盾构开挖过程中，土仓的排土对刀盘前方土体应力稳定有着至关重要的作用.
图5 盾构始发掌子面前方1D截面水平位移曲线
Fig.5 Curves of horizontal displacements apart away 1D from cutter face while shield launching
图6 接触力链图
Fig.6 Force chains distribution
2.3 掌子面应力分布及平均应力
在盾构掌子面位置设置系列测量圆测量掌子面土体的水平应力，当盾构进尺为D/6、D/3、D/2位置时，根据每个测量圆的位置及水平应力值，计算掌子面平均水平应力值，绘制平均水平应力随盾构推进的变化曲线.以掌子面中心为坐标点(0,0)绘出掌子面土体的水平应力分布曲线.Chen等[8]通过离散元研究了盾构隧道的掌子面支护应力问题，与文献[8]相似，绘制的掌子面土体的水平应力分布曲线如图7所示.从图7可以看出，在掘进过程中，掌子面的平均水平应力在不排土与排土情况下，都经历了一个先迅速降低，到一个最低值后又提升至较稳定水平应力的过程.不排土工况下，掌子面位置土压力上升明显，掌子面平均应力在盾构推进约6 cm(D/5)位置开始稳定，基本保持在70 kPa左右；排土工况下，掌子面平均应力在盾构推进约5 cm(D/6)位置开始稳定，基本维持在15 kPa左右.
图7 掌子面平均水平应力变化曲线
Fig.7 Changing curve of average horizontal stress at the tunnel face
如图8所示，不排土工况下，当盾构掘进分别为D/6、D/3、D/2位置时，掌子面平均应力的分布大体相似，掌子面应力值随位置变化的趋势不明显.在排土工况下，土仓排土使得靠近螺机位置掌子面应力非常小；当远离螺机位置时，掌子面应力逐渐增大，在距离螺机一定距离(约D/3位置)，应力开始保持均匀，这说明实践中将刀盘推力除以刀盘面积得到掌子面平均应力，并把此平均应力当作掌子面应力的方法是值得考量的.当挤土效应明显时，此方法可行，但在排土情况下，用此方法得到的掌子面应力存在较大误差.
图8 掌子面上水平应力分布曲线
Fig.8 Distribution curve of horizontal stress at the tunnel face
(1)在不排土工况下，刀盘对盾构上方、前方和下方均有明显的挤密效应，其中以前方挤密最为严重.此时竖向和水平位移在刀盘附近达到最大，向地表和刀盘前方逐渐扩散变小；在排土工况下，刀盘周围土体无明显挤密，水平位移很小，竖向位移在刀盘处较小，向地表延伸逐渐变大，呈现出明显的沉降叠加效应.
(2)在不排土工况下，刀盘前方因挤土产生的竖向位移呈椭圆形，水平位移呈圆形向距刀盘较远处扩散变小；在排土工况下，整个土层水平位移很小，竖向位移呈不规则漏斗状向地表延伸变大.
(3)盾构土仓的排土状态对掌子面前方土体的应力影响较大.在不排土工况下，掌子面前方应力增大明显，在排土工况下，掌子面前方土体应力基本与初始地层一致，这说明地层应力受盾构掘进影响较小.
(4)在不排土与排土工况下，掌子面平均应力随盾构掘进都经历了一个先降低而后升高至稳定值的过程.在不排土工况下，土体水平应力在掌子面的分布基本一致；在排土工况下，由于螺机的存在，土体水平应力分布较为不均，螺机附近土压小，盾顶附近土压大.对于螺机排土情况，实践上以掌子面平均应力来代替整个掌子面应力的方法有待考量.
参考文献：
[1] 魏纲. 盾构法隧道施工引起的土体变形预测[J].岩石力学与工程学报,2009,28(2):418-424.
[2] 梁荣柱,夏唐代,林存刚,等. 盾构推进引起地表变形及深层土体水平位移分析[J]. 岩石力学与工程学报,2015,34(3):583-593.
[3] LI P, DU S J, MA X F, et al. Centrifuge investigation into the effect of new shield tunneling on an existing underlying large-diameter tunnel [J].Tunn Undergr Sp Tech，2014,32:59-66.
[4] 徐前卫,朱合华,廖少明. 砂土地层盾构法施工的模型试验设计研究[J].地下空间与工程学报, 2006, 2(3):361-368.
[5] 蒋洪胜,侯学渊.盾构掘进对隧道周围土层扰动的理论与实测分析[J].岩石力学与工程学报,2003,22(9):1514-1520.
[6] 姜忻良,崔奕,李圆.天津地铁盾构施工地层变形实测及动态模拟[J].岩土力学,2005,26(10):1612-1616.
[7] MAYNAR M, RODRIGUEZ M. Discrete numerical model for analysis of earth pressure balance tunnel excavation[J]. Geotech Geo Eng，2005，131(10):1234-1242.
[8] 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].Comput geotech,2011,38(2):187-195.
[9] JIANG M J, YIN Z Y. Influence of soil conditioning on ground deformation during longitudinal tunneling[J].Comptes rendus mecanique， 2014,342:189-197.
[10]JIANG M J, YIN Z Y. Analysis of stress redistribution in soil and earth pressure on tunnel lining by discrete element method[J].Tunn Undergr Sp Tech,2012,32:251-259.
[11]蒋明镜,王富周,朱合华. 考虑尾隙的盾构隧道土压力离散元数值分析[J].地下空间与工程学报,2010,6(1)：29-32.
[12]王俊,何川,封坤,等. 砂卵石地层中大断面泥水盾构泥膜形态研究[J].现代隧道技术,2014,51(6)：108-115.
[13]江英超,方勇,何川,等.砂卵石地层盾构施工滞后沉降形成的细观研究[J].地下空间与工程学报,2015,11(1)：171-177.
[14]王明年,魏龙海. 卵石地层中地下铁道施工力学的颗粒离散元法模拟技术及应用[M].成都:西南交通大学出版社,2010：33-34.
Abstract：In order to clarify the disturbance mechanism of the soil around the earth pressure balance(EPB) shield affected by the squeezed and exhausted soil situation of tunneling, discrete element method (DEM) was adopted to establish two groups of two-dimensional numerical models of EPB shield tunneling. The comparisons were carried out to analyze the mesomechanical characteristics of soil at the tunnel face and its surrounding strata during the shield tunneling. This research shows that the shield squeezed the surrounding strata obviously when the soil in the chamber was not exhausted. The vertical displacement developed ovally, while the horizontal displacement developed circularly, all of which decreased with an increase in distance along the horizontal direction. When the soil in the chamber was exhausted, the soil ahead of the cutterhead was not densified obviously, the horizontal displacement was low, and the vertical displacement developed to be irregularly funnel-shaped and higher. When the soil in chamber was not exhausted, the influence of shield tunneling on soil stress ahead of the tunnel face was significant, and the distribution of the average horizontal stress at the tunnel face was uniform. However, with the soil being exhausted, the stress of soil ahead of tunnel face was similar to the initial soil condition, and the distribution of the average horizontal stress at the tunnel face varied greatly due to the soil exhausting.
Key words：dry sand; earth pressure balance shield; soil exhausting; ground response; discrete element method
文章编号：1671-6833(2017)01-0016-06
收稿日期：2015-09-21；
修订日期：2015-11-24
基金项目：国家自然科学基金青年资助项目(51208516)；中南大学研究生创新项目(2016zzts414)
中图分类号：U455.4
文献标志码：A
doi：10.13705/j.issn.1671-6833.2016.04.028