1、引言
早期的隧道开挖模拟常假设为二维平面应变问题，为了反映隧道开挖面引起的应力松弛，学者们提出了不同的处理方式，如convergence-confinement method（或Stress relief method）（图1），progressive softening method（图2）, gap method（图3）等。这些方法没有完整体现开挖的三维效应。
图1
图2
图3
隧道开挖的三维模拟一般通过“step-by-step”的形式实现，即首先无支撑下开挖一段土体（单元移除），然后在相应位置设置衬砌管片，然后进行下一步开挖。近年来也有学者尝试完整模拟真实的开挖过程，如开挖机具（盾构）推进、灌浆等https://doi.org/10.1002/nag.395)。本系列首先介绍step-by-step分步三维开挖模拟方法。
2、算例要点
（1）考虑弹性模量随深度变化（Gibson地基）：利用材料参数的温度相关性。
（2）利用蒙皮（skin）模拟衬砌。
（3）利用单元生死模拟隧道施工。
3、算例概况
该算例改编于Franzius J N. Behaviour of buildings due to tunnel inducedsubsidence[D]. Department of Civil and Enviromental Engineering. ImperialCollege of Science, Technology and Medicine. London, SW7 2BU, 2003.利用对称性，模型如图4所示，在地面以下30.5m开挖一个直径4.75m的隧道，开挖沿x负向进行。隧道长度100m，模型x向总长155m，y向总长80m，z向总长80m。每段开挖2.5m，共40步挖完，本例只给出前20步的结果。
土体采用莫尔库伦模型，凝聚力5kPa，摩擦角25度，剪胀角12.5度，地下水位于地表，有效重度10kN/m3，水平土压力系数1.0，弹性模量随深度向下线性增加，即E=4000（1+0.2*depth）。
衬砌混凝土采用弹性模型模拟，参数如图5所示。
图4
图5
4、建模过程
（1）通过拉伸建立三维实体，partition出隧道位置（开挖土体位置），为划分网格方便，可在隧道周围分割出一个矩形，便于网格过渡（图6）。为便于后期操作，将每步开挖的土体创建集合r1、r2等。
图6
（2）在Property模块中创建土体材料，注意选择参数与温度相关（图7）。实现弹性模量的分布也可以通过参数随场变量变化的功能实现，但是场变量的定义不能在cae中进行。方便起见，这里利用温度分布（可以在cae中设置）来实现。图7中的温度0对应模型底面，60对应模型顶面。建立相应section并分配给土体。
图7
（3）创建衬砌混凝土材料，并创建板shell截面属性（图8），输入厚度0.168。
图8
（4）创建蒙皮。在模型中利用createdisplay group（按钮或tools菜单下）将开挖土体移除，执行special-skin-create将第一段开挖土体后余下土体的表面创建蒙皮（skin）（图9）。
图9
（5）执行Assignsection，将板截面属性分配给蒙皮。注意将选择对象改为skin，在屏幕底端顺便建立集合s1（图10）。依次建立所有的衬砌集合。将所有衬砌建为集合sall。
图10
（6）拼装实体后，进入step模块，建立geostatic分析步，用于初始应力平衡（时间步长可选自动）。然后依次建立r1、s1、r2、s2……等分析步（依次为第一段开挖、第一段支撑、第二段开挖……）。将输出控制中增加temp温度作为输出变量（output菜单）。
（7）在Interaction模块中创建model change，将所有的衬砌蒙皮集合在geo分析步中移除；将r1在r1分析步中移除；将s1在s1分析步中重新激活（图11），余同。若开挖步数很多，在inp文件中直接写语句比较方便。
图11
（8）在load模块中，将模型底部所有方向位移约束住，限值周边相应的水平方向位移。在geo分析步中施加z向体力-10。设置相应的初始应力场（图12）和相应的温度场（图13）。图13中的analytical field-1为定义的分布，可以在tools菜单下定义，也可以直接点击f（x）按钮定义，本例中直接为Z，即温度等于Z坐标，模型顶面温度为60，对应图7中的模量4MPa。
图12
图13
（9）划分网格后提交计算，蒙皮用板单元划分。
5、结果简单分析
首先检查温度设置是否准确（图14）。
图14
图15是第一步开挖后隧道附近土体的变形，隧道底部上抬，隧道顶部下沉。在平面分析中，设置了衬砌后的开挖引起的底部隆起会通过刚性衬砌上传到顶部，造成沉降偏小，甚至有可能出现土体上抬变形的情况。三维分析可以改善这一情况。考虑随深度增加的弹性模量也会缓解上抬变形。
图15
图16是x=0，z=60（地表）开挖过程中地表竖向位移的分布，图17是y=0，z=60（地表）开挖过程中地表竖向位移的分布，可以观察到明显的三维效应。图18是地表位移动态变化过程，该图通过overplot功能绘制（图19）。
图16
图17
图18
图19
来源： ABAQUS在岩土工程中的应用
微信
QQ
微博