FLAC 3D软件在地铁隧道地震动力响应的应用

张蒙蒙 雷震宇

（同济大学铁道与城市轨道交通研究院，201804，上海∥第一作者，硕士研究生）

1 FLAC 3D软件

FLAC（Fast Lagrangian Analysisof Continua）是由美国ltasca consulting Group Inc.开发的二维显式有限差分法程序，采用适合于模拟大变形问题的拉格朗日方法，即模型单元网格随着变形而不断更新。拉格朗日法按照时间步采用动力松弛的方法来求解运动方程，采用混合离散化方法有效模拟模型材料的塑性破坏和流动。

区别于其他有限元数值模拟方法，FLAC 3D软件使用显式有限差分方法求解。对于显式法，非线性本构关系与线性本构关系没有算法差别，可以方便地求出应力增量、不平衡力，并跟踪系统的演化过程。另外，在有限差分法中，基本方程组和边界条件近似地改为差分方程表示。由于使用显式方法按时间递步解算代数方程避免了有限元方法中隐式矩阵解算微分方程，不需要建立刚度矩阵，省去了求解大型联立方程组的步骤，所以占用计算机内存少，节约了求解时间，提高了计算效率。

FLAC 3D软件采用全动力学运动方程，可以分析和计算非稳定过程，即使在静态模拟中，也利用动态运动方程计算，使得数值模拟不存在力的不稳定过程，消除了数值上的障碍，这是有限元方法不能解决的。

FLAC 3D软件中单元材料可采用线性或非线性本构模型。在外力作用下，当材料发生屈服流动后，网格能相应发生形变及移动，因此，在模拟大变形问题及其他塑性问题时，FLAC 3D软件所采用的混合离散方法比有限元法中通常采用的离散集成法更为准确合理。在动力模块中，FLAC 3D软件采用完全非线性分析方法，适用于任何指定的非线性本构模型；非线性材料定律使不同频率的波之间可以自然地干涉和混合；可以自动计算永久变形；采用合理塑性方程，使得塑性应变增量与应力相联系；可以方便进行不同本构模型的比较；可以同时模拟压缩波和剪切波的传播及两者耦合作用时对材料的影响。

FLAC 3D软件中内置了11种材料本构模型，包括空单元模型、3种弹性模型、7种塑性模型。每种模型对应1种特殊类型岩土材料的本构特性。模型的任何一个个体单元都可以采取不同的本构模型或赋予不同的材料参数。

2 建模

2.1 模型建立

本文依照上海轨道交通8号线青云路段地质状况建模，该段主要的土层分布和土体物理参数如表1所示。地铁隧道外径为6m，隧道衬砌厚度为0.3m；C55混凝土，密度为2 600kg／m3，弹性模量E＝35 500 MPa；埋深分别采用9m、15m、30m。由于FLAC 3D软件采用显式有限差分方法进行数值模拟，相比较于ANSYS等有限元软件，其对于网格尺寸精细度的要求更低，因此可以采用较大的单元尺寸。

表1 各土层岩体物理参数

由文献［1］知，模型尺寸取隧道直径的7～8倍时边界影响已较小，故建模为：模型宽度100m，厚度方向60m，纵向沿隧道方向延伸10m，岩土体采用摩尔－库伦本构模型，隧道衬砌采用弹性体本构模型。另外，根据FLAC 3D软件动力计算中网格的尺寸取决于输入波形的最高频率，通常取最高频率对应波长的1／10～1／8的理论，以及本文所采用的震动波在各岩层波速及波长推算出的模型最大单元尺寸，为方便建模采用1m×1m网格计算，共生成71400个单元，79 948个节点。根据以往的文献和计算经验，这个精度是满足要求的［3］。

2.2 地震波的矫正和输入

FLAC 3D软件输入的动力荷载可以加载在模型边界或内部节点。加载荷载的方式为向模型输入应力时程、加速度时程、速度时程、集中力时程。荷载的输入方式可为FLAC 3D软件内嵌函数Fish语言定义的规律波形，也可以由table表类读入离散波形。

在FLAC 3D软件动力计算中，网格的尺寸取决于输入波形的最高频率，通常取最高频率对应波长λ的1／10～1／8。即输入波的最大频率越高，网格的尺寸要求越小。为了保证地震波在土体中的传播精度，同时控制网格尺寸及数量以节约计算时间，在波形输入之前需要进行滤波，把波形中的高频成分滤掉。另外，对于各个完整波形，输入加速度时程积分得到的最终速度和最终面积均不为0，故动力计算过程结束时会在模型底部出现继续的速度和残余位移，为此需要进行基线校正。SeismoSignal软件是一款利用惠更斯原理积分的波处理软件，可以对加速度时程的1次积分和2次积分结果进行修正。

本次数值模拟为1995年日本阪神地震波，利用SeismoSignal软件将地震波中高于5Hz的波形成分过滤掉，并对加速度时程进行了基线校正。

本文采用从模型底部水平（X）方向和竖直（Z）方向同时加载。经滤波和校准的水平向、竖直向地震波如图1和图2所示。水平向地震波峰值为0.6 g，竖直向地震波峰值为0.2 g，作用时间为20s。

图1 经过滤波和校准的水平向地震波

2.3 确定边界条件

在动力计算中，模型边界会存在波的反射，为了更接近真实的情况和减小数值计算负担，采用合理的模型尺寸和利用合理边界条件吸收入射波是十分重要的。FLAC 3D软件动力分析的边界条件主要分为静力边界（粘性边界）条件和自由场边界条件。对于地震动力荷载，在刚性地基上可采用自由场边界，模型底部不必加静态边界；在柔性地基上需在模型周围施加自由场边界，同时在模型底部施加静态边界条件。

图2 经过滤波和校准的竖直向地震波

本文模型采用自由场边界。自由场边界利用apply ff，apply ff＿corner，apply ff＿side命令在建立的基本模型四周及4个角生成平面网格和柱网格，通过在边界设置切向和法向的独立阻尼器与主体网格进行耦合。自由场边界中的独立阻尼器可有效吸收来自模型内部切向和法向的入射波，模拟真实的无限边界中波和能量的传播效应。切向粘滞力fs和法向粘滞力fn分别为：

式中：

υn——边界上的法向速度分量；

υs——边界上的切向速度分量；

ρ——介质的密度；

CP——P波的波速；

CS——S波的波速。

2.4 确定阻尼

FLAC 3D软件采用动力方程求解准静力问题和动力问题，在计算中设定阻尼是必要的。对于动力问题中的阻尼，需要在数值模拟中重现自然系统在动荷载加载时的阻尼大小。FLAC 3D软件提供了瑞利阻尼（rayleigh damping）、局部阻尼（local damping）和滞后阻尼（hysteretic damping）。本文采用局部阻尼。

局部阻尼在振动循环中通过在节点或结构单元节点上增加或减小质量的方法达到收敛。由于增加的单元质量和减小的单元质量相等，就一次总体来说，系统保持质量守恒。由于局部阻尼系数不需求解系统的自振频率，且相对于瑞利阻尼不会减小时间步，所以在求解小型模型的简单问题中有很大的优势。当节点速度的符号改变时增加节点质量，当速度达到最大值（或最小值）时减小节点质量，因此损失的能量ΔW／W 是最大瞬时应变能Wd的一定比例，ΔW／W 是与频率无关的。因为ΔW／W 是临界阻尼比D的函数：

式中：

αL——局部阻尼系数；

D——临界阻尼比，取值可参考瑞利阻尼的ζmin。

和祥轩的老板以抢占民女罪被县衙抓走。最后的结局是不得不使上一笔不菲的赎金把人给赎回来。和祥轩老板仿佛是死过了一回，他亲眼见证了那沾满血迹的刑具是如何的寒威逼人。他只在那刑房里待了不到一个时辰，竟然就把屎尿全拉在了裤裆里。所以，以后，和祥轩的人只要瞧见那些个穿半头鞋子的人都敬而远之。更要命的是，那案子还久拖而不得了结，和祥轩还得定期给官府送银票续保。

对于岩土材料，临界阻尼比一般取0.05，则局部阻尼系数可取0.157。

3 计算分析

3.1 初始地应力

在施加动力荷载前，需要对岩土层的密度、体积模量、剪切模量、黏聚力、摩擦角等物理参数赋值，并在竖直方向施加重力，solve软件会自动求解原始的地应力情况以还原真实的岩土力学情况。

为了模拟真实的地铁隧道周围应力状况，分步对左右线隧道进行开挖和支护。考虑到地铁隧道管片宽度为1.3m，因此采用的开挖模拟方法为：每次沿隧道轴向利用Null模型向前掘进1.3m，通过Solve软件得到开挖后应力状态，然后对前一环开挖段添加初期支护并再次通过Solve软件计算，直到计算平衡，再逐次向前推进；先对左线开挖、支护，之后依照相同方法对右线开挖并支护，以此得到支护完成后的地应力场（见图3、图4）。

图3 埋深9m时左线开挖支护后竖直应力云图

图4 埋深9m时左右线开挖支护后竖直应力云图

3.2 动力计算和分析

从模型底部由X方向和Z方向在模型X边界和Z边界同时加载处理过的地震波，动力作用时间10 s，记录每1s的最大主应力和最小主应力、隧道及周围土体单元X向和Z向应力、关键点位移和速度。

3.2.1 关键点位移

取左线隧道洞顶衬砌内侧和外侧点、隧道肩部点、隧道左右两侧衬砌内外侧点作为主要监测点（见图5），记录埋深为9m、12m、30m时，隧道顶、肩和两侧监测点在震动输入的竖直方向和水平方向每秒钟的位移情况，如图6～11所示。其中，GP1为隧道顶部衬砌内侧一点；GP2为隧道顶部衬砌外侧围岩上一点；GP5为隧道肩部衬砌外围岩上一点；GP10为隧道衬砌壁远离相邻隧道方向的侧面上一点；GP11为隧道衬砌壁相近相邻隧道方向的侧面上一点。

图5 左线隧道各考察点位置

图6 埋深9m时各考察点竖直位移变化

图7 埋深9m时各考察点水平位移变化

图8 埋深12m时各监测点竖直位移变化

图9 埋深12m时各监测点水平位移变化

图10 埋深30m时各监测点竖直位移变化

图11 埋深30m时各监测点水平位移变化

从图中可以看出，埋深不同，对于各点在竖直和水平地震波传播方向的位移影响较大，随着埋深增加，隧道顶、肩、两侧位移更加趋于一致，位移量减小，隧道洞壁和衬砌变形变小，更加不易破坏。这说明：随着上覆土层体积的增加，重力及周围岩体性质对隧道衬砌的约束增强，使得隧道衬砌对地震动力的响应减弱。

另外，隧道衬砌在整个地震动力响应过程中的最大位移响应出现在第5s，这为之后的分析提供了依据。

图12～17表示在位移响应最大的第5s时隧道在埋深分别为9m、12m、30m情况下的竖直应力云图和水平应力云图（单位：Pa）。

表2是隧道顶、肩、两侧考察点周围单元在位移值峰值5s时，埋深分别为9m、12m、30m情况下的竖直应力值（szz）和水平应力值（sxx），以及埋深增加后（9m→12m、12m→30m）应力值增长百分比。

图12 埋深9m，5s时竖直应力云图

图13 埋深9m，5s时水平应力云图

图14 埋深12m，5s时竖直应力云图

图15 埋深12m，5s时水平应力云图

图16 埋深30m，5s时竖直应力云图

图17 埋深30m，5s时水平应力云图

由图12～17及表2可以看出：① 在埋深增加的情况下，隧道顶部在竖直方向拉应力增长，压应力区向两肩扩散，且数值逐渐增大；② 埋深越大，隧道单元应力值越大，但增长百分比越小，增长幅度明显减小。这说明：在埋深较深的地方围岩对隧道的约束越大，地震波对隧道衬砌的内力影响相对越弱。

3.2.3 变形云图

FLAC 3D软件的后处理功能可以绘出模型的变形云图，用以形象地表征变形的发生方向和相对大小。图18～20采用模型的变形云图、变形后网格叠加的方式显示了埋深分别为9m、12m、30m，震动持 续5s时（位移变形最大时刻）的位移变形情况。

表2 不同埋深下隧道各考察点周围单元应力值

图18 埋深9m变形云图

图19 埋深12m变形云图

图20 埋深30m变形云图

从图中可以看出，上覆土层厚度为12m时的隧道变形情况较埋深为9m、30m时更加严重。这主要是由于隧道所处土层造成的，此时隧道处于淤泥质粉砂黏土层和黏土层，相对于其他情况所处的粉质黏土层，其弹性模量E、摩擦角f及内聚力c较小，由于周围岩土控制地震波的传递和衰减特性，地震波对于软弱地层造成的影响显然比坚硬地层大得多。

4 结语

本文通过FLAC 3D软件模拟地铁双线隧道在不同埋深情况下，在竖直方向和水平方向同时加载地震波时的应力和位移变形情况，得出以下结论：

（1）上覆土层厚度12m时，隧道的变形情况比其他情况下更加严重。由于隧道处于淤泥质粉砂黏土层和黏土层，其E，f，c相对于其他情况所处粉质黏土层较小，根据周围岩土控制地震波的传递和衰减的特性，因此软弱地层地震波造成的影响显然要大的多。

（2）在埋深增加的情况下，隧道顶部在竖直方向拉应力增长，压应力区向两肩扩散，且数值逐渐增大。

（3）埋深越大，隧道单元应力值越大，但增长百分比越小，增长幅度明显减小，说明在埋深较深的地方围岩对隧道的约束越大，地震波对隧道衬砌的内力影响越弱。

（4）埋深不同，对于各点在竖直和水平地震波传播方向的位移影响较大，随着埋深增加，隧道顶、肩、两侧位移更加趋于一致，位移量减小，隧道洞壁和衬砌变形变小，更加不易破坏。

综合各结果可知，随着隧道上覆土层厚度的增加，隧道的破坏程度逐渐减小。这是由于土体重度增加和对隧道衬砌产生影响的岩体量增加，约束增强，使得地震波的动力作用对隧道衬砌的位移和内力影响减弱。

利用FLAC 3D软件对地铁隧道进行初始应力场、开挖和简单支护、地震波作用模拟，可直观地体现地铁隧道在动力作用下围岩位移和应力分布情况，并据此进行一系列探讨。相比其它有限元分析软件，FLAC 3D软件更加节约计算机内存和计算时间，其强大的后处理功能使分析更加直观和易得。因此，FLAC 3D软件在隧道地震动力响应方面的应用具有显而易见的适用性和优越性。由于时间等因素的限制，本文所采取的模型仍有许多值得改进之处：① 模型底部边界若采用静态边界限制或可得到更精确的分析结果；② 隧道的衬砌只采用简单的初期支护，或可采用更加精细的支护模拟以更加接近真实工况；③ 本文采取震动每秒结束时的位移和应力情况进行分析，或可采用峰值进行分析对比。这些问题将在以后的研究工作中逐步完善。

［1］王铁男.FLAC在地铁隧道数值模拟中的应用［J］.沈阳大学学报，2010，22（1）：11.

［2］唐益群，栾长青.地铁振动荷载作用下隧道土体变形数值模拟［J］.地下空间与工程学报，2008，4（1）：105.

［3］陈育民，徐鼎平.FLAC／FLAC 3D基础与工程实例［M］.北京：中国水利水电出版社，2011：75－88，297－312.

［4］彭文斌.FLAC 3D实用教程［M］.北京：机械工业出版社，2011：207－249.

［5］李围.隧道及地下工程FLAC解析方法［M］.北京：中国水利水电出版社，2009：238－245.

［6］陈贵红.沉管隧道抗震数值分析［D］.成都：西南交通大学，2002.

［7］邬玉斌，刘如山，杨学山.地铁交通队某实验楼的环境振动的模拟分析［J］.世界地震工程，2008，24（4）：70.

［8］马乐元，丁毅，郭传创.基于FLAC 3D的移动荷载作用下基坑结构响应的分析［J］.工程与建设，2011，25（1）：90.

［9］曹恒将.考虑动载的巷道顶板离层及控制的FLAC 3D模拟研究［D］.青岛：山东科技大学，2010.

［10］刘波，叶圣国，陶龙光.地铁盾构施工引起邻近基础沉降的FLAC元数值模拟［J］.煤炭科学技术，2002，30（10）：9.

［11］苏丽娟.高速铁路隧道围岩支护参数优化设计［D］.北京：北京交通大学，2011.