张义同,袁 杰
(天津大学 机械工程学院力学系,天津300072)
在地铁盾构施工中,施工事故时有发生,其中社会影响较大的事故是由掘进面失稳引发地表沉降导致地表坍塌造成的[1]。影响掘进面稳定性的因素很多,土体的应变软化可能引起隧道掘进面失稳,但在隧道开挖面稳定性分析中还很少考虑应变软化的影响。早在20世纪60年代,人们就开始关注土体的应变软化行为,Skempton在研究边坡的稳定性时就指出了,具有应变软化性质土坡与具有应变硬化性质土坡的渐进破坏过程不同,应变软化对边坡的稳定具有重要影响[2-3];张嘎基于瑞典条分法的基本原理,提出了一个能考虑土的抗剪能力与剪应变关系的应变软化边坡稳定性简化评价方法,研究了土的应变软化特性对边坡稳定性的影响[4]。Brown较早地考虑了涉及应变软化时隧道围岩稳定性分析的理论问题,但没有考虑塑性区弹性应变的变化,得到的结果存在较大误差[5];Lee和Pietruszczak对隧道围岩的应变软化进行分析,采用了有限差分法求解平衡方程与协调方程,得到了隧道围岩的塑性区与弹性区的分布范围[6]。于学馥利用材料的参数从弹塑性界面向隧道壁减小来考虑材料软化,分析了隧道围岩的稳定性[7];王星华则采用弹性 塑性软化 塑性残余三线性应力 应变模型推导了隧道围岩的弹塑性区半径、围岩应力的解析表达式,分析表明应变软化对围岩稳定产生不利的影响[8]。
研究表明,应变局部化是应变软化材料的一种失稳现象[9]。应变局部化的变形往往集中在一个很狭窄的带状区域内,也称为剪切带。目前已有一些文献对土体剪切带进行了研究[10-14],在理论研究中,一般把应变局部化看成是一个分岔问题[14-15],Zhang[9]进一步指出,其稳定的分岔解是一个两相解,用一般的分岔分析方法求出的是单相解,难以求出其稳定的分岔解,并通过相变分析,给出了土体平面应变下剪切带的解析解。盾构掘进面的失稳是三维问题,对其进行解析分析是比较困难的,采用数值模拟方法对盾构掘进面的稳定性进行分析,对比土体是否考虑应变软化的情况,得出应变软化对掘进面稳定产生的影响。
首先介绍利用相变理论对Alshibli和Sture[10]实验中砂土剪切带的解析分析,再利用数值方法对剪切带的形成与扩展进行数值模拟,将模拟结果分别与解析解和实验结果进行比较,来检验本文数值方法的正确性。
Abeyaratne[16]和 Zhang等[9]以具有应变软化行为的一维弹塑性杆为例,说明了在均匀拉伸的过程中可能存在的两个分岔解:一个是传统的单相均匀变形解,另一个是两相解。并指出两相解的应变能最小,是稳定的分岔解。Fu[17]和Zhang等[9]建立了应变局部化的两相平衡模型,剪切带内区域是高应变相(塑性相),剪切带以外的区域是低应变相(弹性相),在每一相内,必须满足平衡方程(不考虑体积力):
式中:div为散度算子;π为第一Piola-Kirchhoff应力张量;T表示张量的转置。
相变界面处的Eshelby力和位移是连续的,应变梯度和应力是不连续的,跨越相变界面的跳跃应满足:
土体剪切带的变形是非均匀变形,可以将剪切带内和剪切带外的变形看作是分片均匀变形,于是各相域内的平衡方程自然满足,由(2)、(3)式可得到平面应变状态下砂土剪切带的控制方程,解控制方程与=1组成的方程组可求出剪切带出现时的临界应力、剪切带的法线方向矢量、剪切带内和剪切带外的变形等。Zhang等[9]对 Alshibli和Sture实验中F4砂土剪切带的解析分析的临界应力和剪切带倾角列于表2中。
采用有限差分法程序FLAC3D进行数值分析,本构模型取为Mohr-Coulomb应变软化模型,材料塑性屈服后,粘聚力、摩擦角、剪胀角发生应变软化,把这些变量作为塑性应变的函数。试件使用Alshibli和Sture实验中的F4砂土,其材料参数如表1所示。
表1 F4砂土材料参数
砂土试样宽0.15m,高0.3m,厚0.02m,共划分900个单元,为了激发局部变形的产生,在试样中引入缺陷单元。边界条件如下:在试样上下两端施加常速率1.5×10-6m/时步以控制位移,左右两侧施加15kPa围压,前后两侧约束位移。数值模拟出砂土试样的变形情况如图1所示,变形从缺陷点启动,缺陷点及其相邻单元中的应力与应变关系表现为应变软化,而远离缺陷点的单元还处于弹性状态,变形充分发展后,集中于局部区域并沿局部区域迅速扩展,形成了一条明显的剪切带(图1(c)),图2为Alshibli和Sture实验的剪切带照片。通过查网格法计算出图1中剪切带的切向与水平方向的夹角为57.9°,剪切带形成时的临界应力为177kPa,与解析分析的临界应力172.056kPa比较接近,也在上屈服极限239.25kPa之下,与实验观测的结果一致,在上屈服极限就要到达之前,剪切带就出现了。剪切带的实验结果、解析解和数值模拟结果对比情况如表2所示,数值模拟的剪切带倾与实验测得的剪切带倾角相吻合,误差仅为5.3%,与解析分析的剪切带倾角误差也仅为4.8%,说明了用本文的数值方法来模拟应变软化土体的变形是可行的。
图1 剪切应变增量发展云图
图2 实验中剪切带的照片
表2 剪切带的数值模拟结果与实验结果和解析解的对比
开挖隧道是一个逐渐推进的过程,考虑到分析的重点是土体应变软化对掘进面稳定性的影响,因此忽略注浆压力、注浆厚度等因素的影响。数值模拟时一次性开挖隧道一定距离,并及时设置衬砌管片,同时在开挖面上施加大小与原始地层侧向静止土压力值相等的梯形支护作用力,地层水平侧向土压力为垂直压力乘以侧压力系数,侧向压力系数取值0.5。为了方便描述,取隧道中心点支护压力值来代表开挖面支护压力大小,支护压力开挖面中心点施加的支护压力,p0为中心点原始地层的水平静止土压力之比。
隧道计算模型如图3所示。模型总长度100m,宽度50m,总高度65m,隧道直径为10m,隧道埋深20m,隧道已开挖50m。模型的坐标如下:y轴沿隧道轴线,z轴垂直于地面。考虑到土体模型对称性,对其一半进行网格剖分。其边界约束条件为:1)地表面为自由面;2)四周限制水平位移,底部限制竖向位移。假定地下水位与地表一致,土体材料为砂土,考虑砂土应变软化时采用Mohr-Coulomb应变软化模型,不考虑砂土应变软化时采用Mohr-Coulomb模型,砂土参数如表3所示,管片材料为C50弹性钢筋混凝土结构,厚度为35cm,采用shell单元模拟。
图3 隧道计算模型
表3 土体材料参数
2.2.1 不考虑渗流时土体应变软化与不软化的对比 当开挖面支护压力小于地层原始地应力时,引起土体应力释放,导致掘进面土体变形。不考虑渗流时掘进面上中心点的水平位移与支护压力比的关系如图4(a)所示,随着支护压力比不断减小,开挖面中心点的水平位移量逐渐增大,当支护压力比减小到一定程度时,中心点的水平位移量有一个突变,之后位移急剧增大,在支护应力几乎不变的情况下变形还会继续发展。根据稳定性理论,把开挖面支护压力的微小变化导致开挖面中心点水平位移突变时的支护压力定为极限支护压力,相应的支护压力比称为极限支护压力比。由图可知,考虑砂土应变软化时的极限支护压力比为0.2,大于不考虑软化时的极限支护压力比0.1,相应的考虑软化时的极限支护压力也大于不考虑软化的极限支护压力。
图4 支护压力比和开挖面中心的水平位移关系
2.2.2 考虑渗流时土体应变软化与不软化的对比 由图4(b)可以看出,考虑渗流及土体应变软化时的极限支护压力比为0.45,也大于考虑渗流但不考虑应变软化时的极限支护压力比0.32,土体的应变软化可以使掘进面的极限支护压力变大。与图4(a)相比可知,渗流对掘进面稳定性的影响也比较大,因此,在有地下水存在的情况下不可忽略渗流作用的影响。
支护压力比同取0.45时,渗流情况下Z向(垂直地面方向)的位移云图见图5,不考虑应变软化时,仅在掘进面附近的位移值较大,地表面的Z向位移即地表沉降还不是很大;考虑应变软化时,位移较大的区域已经从掘进面处扩展到地表面附近,在开挖面前方形成一个“漏斗”状,两者的塑性区发展情况见图6。
图5 Z向位移云图
图6 土体的塑性区云图
从图6可以看出,两者的塑性区扩展范围差别很大,不考虑砂土应变软化时塑性区虽然得以扩展,但扩展范围还比较小,还没有发展到地表面;而考虑应变软化时,掘进面前方土体的塑性区已迅速大范围扩展,并与地面贯通,引起地表产生较大的变形。地表沉降情况如图7所示,其中(a)是掘进面正上方地面沿x方向(垂直隧道轴线方向)的沉降曲线,(b)为隧道正上方沿y方向(隧道轴线方向)的地表沉降曲线,开挖方向沿坐标轴正向,开挖面位于横坐标原点处。由图可知,考虑到土体的软化时,隧道正上方地表已经下沉140mm见图7(a),地表的最大沉降出现在开挖面前方,其值为200mm见图7(b),此时地表已塌陷;而不考虑软化时,最大地表沉降仅为12mm。从图还可以看出,土体应变软化对掘进面上方地表附近的沉降影响很大,对远离掘进面的地表影响较小。
图7 地表沉降曲线
由以上对分析可知,土体的应变软化对盾构隧道掘进面的稳定性影响非常大,随着盾构掘进的进行,由于开挖面支撑处理不当,具有应变软化特性的土层会逐渐因局部区域的变形急剧发展而加快盾构掘进面的失稳,导致较大的地面沉降变形,严重时导致隧道坍塌。而且这种局部化变形的发展有时会很突然,没有任何警示的迹象。所以,在隧道工程中,遇到具有应变软化的特性的地质,必须考虑应变软化对掘进面稳定性的影响,从而在施工中采取合理的支护方法,预防隧道土体坍塌事故的发生。
2.2.3 砂土不同应变软化行为对掘进面稳定性的影响 砂土的软化参数剪胀角和内摩擦角取值见表4,图8(a)为不同方案剪胀角的支护压力比与开挖面最大水平位移曲线,随着剪胀角峰值的减小,开挖面极限支护压力比逐渐增大。图8(b)为不同方案内摩擦角的支护压力比与开挖面最大水平位移曲线,随着内摩擦角峰值的减小,极限支护压力比也逐渐增大,相比较而言,内摩擦角的变化对掘进面稳定性的影响较大。
表4 不同方案的剪胀角与内摩擦角
图8 支护压力比与开挖面最大水平位移关系
在天津地铁九号线十一经路站到大直沽西路站区间,采用小松公司的TM634PMX土压平衡盾构机,刀盘直径6.32m,隧道衬砌采用35cm厚的C50混凝土管片,管片环宽1.2m。此区间地质构造复杂,地层复杂多变,其中500~600环的地质剖面图见图9,土层有填土、粘土、粉质粘土、粉土、粉砂。隧道从中部较厚的粉质粘土层穿过,其埋深为14.1~14.5m。
图9 地质剖面图
由于FLAC3D前处理功能较弱,为了建立隧道的有限差分模型,先通过AutoCAD建立三维模型,然后利用ANSYS对模型进行网格划分,最后通过转换程序将模型导入FLAC3D中,从而实现快速建立复杂三维地质模型(图10),模型长120m,宽30m,高度取35.5m,隧道直径6.32m,埋深取平均值14.3m。边界条件为:模型侧面和底面为位移边界,侧面限制水平位移;底部为固定边界,限制水平移动和垂直位移;模型上面为地表面,为自由边界。土体的本构模型采用Mohr-Coulomb应变软化模型,由于有地下水的存在,同时采用水土耦合计算,考虑了渗流作用。各土层的主要物理参数见表5,衬砌管片的密度为2 450kg/m3,弹性模量为34 500MPa,泊松比为0.17。
图10 计算模型
表5 土层的参数
盾构施工引起的地表沉降因素众多,其中开挖面前方土体的沉降或隆起大小与开挖面支护压力的大小密切相关。若支护压力过大,造成掘进面前方地表的隆起;若支护压力过小,引起前方地表沉陷,严重时会引起地表坍塌。图11给出了开挖10m、20m、50m、80m、90m时不同支护压力比下的地表纵向 (沿隧道轴线方向)沉降曲线,开挖方向沿横坐标轴负向。
图11 不同支护压力比下的地表纵向沉降曲线
从图11可以看出,开挖面支护压力比由1.0逐渐减小时将会引起开挖面前方地表一定距离内沉降变形越来越大,最大沉降点发生在开挖面前方5m左右,接近隧道直径6.32m;支护压力比逐渐增大时,开挖面前方地表的隆起量会越大越大,最大隆起点则在开挖面前方大约10m处。开挖10m时,支护压力比为0.7的开挖面前方地表最大沉降值为2.90mm,0.6时最大沉降值达到4.38mm(图11(a))。由于天津地铁九号线要穿越城市中心地带,地面建筑物密集、地下管网密布等,因而施工中对地表的控制要求较为严格,要求盾构掘进时开挖面前方地表沉降不超过3mm,以及隆起量控制在2~3mm。根据要求由图11可得出,开挖10m、20m、50m、80m、90m时掘进面合理支护压力比范围分别为0.6~1.4、0.5~1.4、0.6~1.4、0.7~1.3、0.65~1.3,各开挖面土层原始侧向土压力分别为0.23MPa、0.23MPa、0.22MPa、0.25MPa、0.24MPa,从而得出合理的支护 压 力 范 围 分 别 是 0.092~0.322MPa、0.132~0.308MPa、0.175~0.325MPa。
在盾构施工过程中,土仓压力左右测点分布如图12所示,图13中脉冲式的波动曲线给出了500~600环区间左右测点的土压监测值,并根据上文预测的范围给出了支护压力的上下限,理想支护压力为开挖面原始侧向土压力。监测的土压力值在理想支护压力附近波动,波动的范围与预测的范围基本一致,说明了分析预测的支护压力范围是合理的,考虑土体软化时对隧道掘进面土体变形分析是比较准确的,如果忽略土体的应变软化行为则偏于危险。
图12 支护力分布与测点位置图
图13 土压监测值
通过数值模拟分析,并以工程实例进行了验证,得到以下结果:
1)土体的应变软化可以加快土体塑性区的扩展,使土层及地表的变形加剧,掘进面的极限支护压力也增大,更容易引起掘进面失稳并导致隧道坍塌。
2)不同剪胀角和内摩擦角的应变软化行为对掘进面稳定性的影响会不同,其中内摩擦角的影响更显著。
3)在工程实例分析中,考虑了土体的应变软化和地下水的渗流,预测了施工中的支护压力,预测的范围与实际监测的土仓压力波动范围基本吻合,得到了工程的验证。这为我国目前正大规模进行的隧道工程建设的设计与施工提供依据,对于掘进面的稳定性控制具有重要的现实意义。
[1]刘辉,张智超,王林娟.2004-2008年我国隧道施工事故统计分析[J].中国安全科学学报,2010,20(1):96-100.LIU HUI,ZHANG ZHI-CHAO,WANG LIN-JUAN.Statistical analysis on safety accidents of tunnel construction from 2004 to 2008in China[J].China Safety Science Journal,2010,20(1):96-100.
[2]SKEMPTON A W.Long-term stability of clay slopes [J].Geotechnique,1964,14(1):77-101.
[3]SKEMPTON A W.Residual strength of clay in landslide,fol-ded strata and the laboratory test[J].Geotechnique,1985,35(1):1-18.
[4]张嘎,张建民.基于瑞典条分法的应变软化边坡稳定性评价方法[J].岩土力学,2007,28(1):12-16.ZHANG GA,ZHANG JIAN-MIN.Stability evaluation of strain-softening slope based on Swedish slice method[J].Rock and Soil Mechanics,2007,28(1):12-16.
[5]BROWN E T,BRAY J W,LADANYI B,et al.Ground response curves for rock tunnel[J].Journal of Geotechnical Engineering,1983,109:15-39.
[6]LEE Y K,PIETRUSZCZAK S.A new numerical procedure for elasto-plastic analysis of a circular opening excavated in a strainsoftening rock mass [J].Tunneling and Underground Space Technology,2008,(23):588-599.
[7]于学馥,等.地下工程围岩稳定性分析[M].北京:煤炭工业出版社,1983.
[8]王星华,章敏,王随新.考虑渗流及软化的海底隧道围岩弹塑性分析[J].岩土力学,2009,30(11):3267-3272.WANG XING-HUA,ZHANG MIN,WANG SUI-XIN.Elastoplastic analysis of surrounding rocks of subsea tunnel with consideration of seepage and material softening[J].Rock and Soil Mechanics,2009,30(11):3267-3272.
[9]ZHANG YI-TONG,Qi DE-XUAN,DU RU-XU,et al.Prediction of shear bands in sand based on granular flow model and two-phase equilibrium [J].Journal of Central South University Technology,2008,15(s1):316-321.
[10]ALSHIBLI K A,STURE S.Shear band formation in plane strain experiments of sand [J].Journal of Geotechnical and Geoenvironmental Engineering,2000,126(6):495-503.
[11]HUANG M S,LU X L,QIAN J G.Non-coaxial elasto-plasticity model and bifurcation of shear banding in sands[J].International Journal for Numerical and Analytical Methods in Geomechanics,2010,34(9):906-919.
[12]ROSCOE K H.The influence of strains in soil mechanics[J].Geotechnique,1970,20(2):129-170.
[13]ARTHUR J F R,DUNSTAN T,AL-ANI Q A J,et al.Plastic deformation and failure in granular media[J].Geotechnique,1977,27(1):53-74.
[14]DESRUES J,VIGGIANI G.Strain localization in sand:an overview of the experimental results obtained in Grenoble using stereophotogrammetry[J].International Journal for Numerical and Analytical Methods in Geomechanics,2004,28:279-321.
[15]SULEM J.Bifurcation theory and localization phenomena[J].European Journal of Environmental and Civil Engineering,2010,14(8-9):989-1009.
[16]ABEYARATNE R,BHATTACHARYA K,KNOWELES J K.Stress-work functions with multiple local minima:modeling phase transformation using finite thermo-elasticity[M].∥ FU Y B,OGDEN R W.Non-linear elasticity:theory and applications.Cambridge,UK:Cambridge University Press,2001,431-490.
[17]FU Y B,ZHANG Y T.Continuum-mechanical modeling of kink-band formation in fibre-reinforced composites[J].International Journal of Solids and Structures,2006(43):3306-3323.