多隧道开挖诱发山体变形与应力数值模拟

Numerical Simulation of Mountain Deformation and Stress Induced by Multi Tunnel Excavation

摘要: 为了研究多隧道开挖过程中山体变形规律，借助有限元软件COMSOL Multiphysics建立了多隧道的二维模型，采用弹塑性方法及Drucker-Prager屈服准则。对多隧道穿越山体的情况下进行数值模拟分析，首先研究了隧道开挖前土体的应力状态；然后，进一步研究隧道开挖后土体的弹塑性行为。模拟过程中，将第一步中的土体当作弹性土，在第二步中加入土的塑性材料模型Drucker-Prager。在二维平面应变下求解，得到了隧道开挖应力主要集中在洞口周围，造成的土体变形主要位于隧道上方和左右两侧，变形主要以竖向地表沉降为主，关于隧道轴线对称两侧水平与竖直方向的变形规律相同。模拟结果表明隧道周围塑性区域的地表沉降和宽度，是预测开挖过程中加固材料用量的重要参数。

Abstract: In order to study the deformation law of mountain in the process of multi tunnel excavation, a two-dimensional model of multi-tunnel is established by using the finite element software COMSOL multiphysics. The elastoplastic method and Drucker Prager Yield Criterion are adopted. Firstly, the stress state of soil before tunnel excavation is studied; secondly, the elastoplastic behavior of soil after tunnel excavation is further studied. In the simulation process, the soil in the first step is regarded as elastic soil, and the plastic material model Drucker Prager is added in the second step. The results show that the tunnel excavation stress is mainly concentrated around the tunnel portal, and the soil deformation is mainly located at the top and left and right sides of the tunnel. The deformation is mainly vertical surface settlement. The deformation law of the horizontal and vertical directions on both sides of the tunnel axis symmetry is the same. The simulation results show that the surface settlement and width of the plastic zone around the tunnel are important parameters to predict the amount of reinforcement material during excavation.

1. 引言

目前有限元法已成为隧道开挖最常用的数值方法之一，多数计算软件都是基于有限元法编制的 [1]。在数值模拟方面，杨晓泉 [2] 等利用COMSOL对浅埋段隧道开挖后拱顶和两侧拱腰围岩变形进行了模拟；段军朝 [3] 等利用ABAQUS软件对大跨度双联拱隧道进行了二维数值模拟，并分析隧道所受应力、应变及塑性变形的情况；魏华 [4] 等基于ANSYS对邻近矿洞隧道开挖进行分析，得出隧道开挖后拱顶和拱脚处的应力增大；杨笑天 [5] 等利用MATLAB对隧道开挖后拱顶位移变化进行了模拟；李俊 [6] 通过建立二维有限元模型对联拱隧道的围岩进行稳定性分析，提出不同施工方法对联拱隧道开挖和设计具有一定的指导意义；胡亮 [7] 等模拟浅埋段隧道开挖后施工地层的沉降变形。P. A. Vermeer [8] 等得出施工工艺是准确预测隧道沉降、水平变形和衬砌力的重要因素；Daniel Dias [9] 等采用FLAC^3D对两隧道结构受力和周围地面位移场的发展进行了讨论，结果表明双隧道同时开挖会导致隧道上方沉降较多。以上研究大都是对双拱隧道开挖后的沉降或水平变形单方面的数值模拟研究，未对隧洞开挖前的状况进行较准确模拟，即隧洞开挖后的状况是在开挖前山体实际情况的基础上进行预测分析。

COMSOL可以很好地模拟岩土体的力学性能，还可考虑非线性应力–应变关系，使得实际情况在计算中得到较好的反映 [10] [11]。本文用COMSOL软件进行隧道开挖过程的模拟，对山体开挖前受自重及其他荷载影响后的状态进行了模拟，为模拟隧洞开挖后提供了更符合实际的研究对象。通过对开挖前后不同受力状态的模拟，得出了开挖前后山体的应力分布和应变，以及导致山体产生水平和竖直位移的变化规律。重点对比单隧道和多隧道开挖后山体的位移和应力的变化情况，为隧道的设计与施工提供了参考条件。

2. 有限元法基本理论

2.1. 基本理论

用有限元进行分析时，首先建立山体的几何模型，然后给定边界条件、荷载条件和材料特性。然后采用弹塑性方法及Drucker-Prager屈服准则，求解稳态结果。最后得到应力、应变、位移等结果 [12] [13]。

2.1.1. Mohr-Coulomb准则 [14]

Mohr-Coulomb条件是一种剪应力屈服条件，其具体表达式为：

$f=\frac{I_1\left({\sigma }_{ij}\right)\sin \phi }{3}\left(\cos {\vartheta }_L+\frac{\sin {\vartheta }_L\sin \phi }{\sqrt{3}}\right)\sqrt{I_2\left(S_{ij}\right)}+C\sin \phi =0$ (1)

式中θ_L为应力Lode角：

${\vartheta }_L=\frac{1}{3}{\sin }^{-1}\left(\frac{3\sqrt{3}}{2}\frac{I_3(S_{ij})}{I_2{\left(S_{ij}\right)}^3}\right)$ (2)

Mohr-Coulomb准则在应力空间是屈服面是一个不规则的六棱锥面。

2.1.2. Drucker-Prager准则 [15] [16]

D. C. Drucker和W. Prager在1952年提出一个内切于Mohr-Coulomb六棱锥的圆锥形屈服面，显然它是Mohr-Coulomb屈服条件的下限，则Drucker-Prager准则可表示为：

$f=\alpha I_1\left({\alpha }_{ij}\right)+\sqrt{I_2\left(S_{ij}\right)}+k=0$ (3)

式中 $\alpha$ 、 $k$ 为常数，

$\alpha =\frac{\sqrt{3}\sin \phi }{\sqrt{3}+{\sin }^2\phi }$ (4)

$k=\frac{\sqrt{3}C\cos \phi }{\sqrt{3}+{\sin }^2\phi }$ (5)

3. 有限元仿真分析

此示例模拟隧道开挖期间土壤的行为。地表沉降和隧道周边塑性区宽度是预测隧道开挖过程中需要加固的重要参数 [17] [18]。

Figure 1. Mountain model

图1. 山体模型

为了计算地应力，使用两个研究步骤。在第一项研究中，计算了隧道开挖前土壤的应力状态。在第二个研究中，计算了土壤移除后的弹塑性行为。这要求将第一步计算的应力响应纳入其中。

为了加快计算速度，第一步考虑土的弹性，第二步增加土的塑性材料模型Drucker-Prager。算例采用二维平面应变法求解。以下案例是以中心线为对称轴，取山体的右半部分作为研究对象。

3.1. 山体模型

隧道开挖前山体模型建立。右半部分山体高度为45米，宽度为90米。以山体左下角点为坐标原点(0, 0)。模型如图1所示。

3.2. 单隧道模型

单隧道开挖后的模型。几何结构由45 m深、90 m宽的土层组成。直径为10 m的隧道位于地表以下20 m的对称轴处。地表以下45 m的基岩约束了垂直方向的位移，并用滚柱边界模拟了土壤在水平方向的无限延伸。模型如图2所示。

Figure 2. Dimensions and boundary conditions of a single tunnel excavation example

图2. 单隧道开挖示例的尺寸和边界条件

3.3. 三隧道模型

三隧道开挖后的模型。基于单隧道模型基础上，直径为10 m的隧道位于20 m处，以(45, 25)为圆心。模型如图3所示。

4. 计算条件

计算条件包括模型的材料参数、边界约束、荷载、与屈服准则等内容。

4.1. 材料参数与基本假定

4.1.1. 材料参数

杨氏模量，E = 12 MPa，泊松比v = 0.495。内聚力 = 130 KPa，内摩擦角φ = 30˚，使用Drucker-Prage准则，并将材料参数与Mohr-Coulomb准则相匹配。材料属性见表1。

Figure 3. Dimensions and boundary conditions of three tunnel excavation examples

图3. 三隧道开挖示例的尺寸和边界条件

Table 1. Properties of each material parameter

表1. 各材料参数的属性

4.1.2. 假设条件

由于场地条件的复杂性，对模型作适当简化处理。为了在有限元分析中使问题得以简化又能放映问题的主要特征，在有限元模拟中作如下设定：

①隧道围岩、山体各自均为均质、各向同性体；

②不考虑变形的时效性，即认为所有变形一步到位；

③围岩为弹塑性体，符合广义胡克定律；

④岩土体为二维连续介质；

⑤为了更好的分析隧道的开挖效应，计算过程中不考虑支护。

4.2. 边界约束与荷载

· 在下边界用固定的约束边界固定位移。

· 在左边界使用对称，在右边界使用滚轮。

· 在顶部保留默认的自由边界，从该区域清除土壤后，在隧道壁上保持自由边界。

· 添加重力节点以考虑重力效果。

4.3. 边界约束与荷载

采用Drucker-Prager (德鲁克–普拉格)准则。

5. 模拟结果与分析

5.1. 模拟结果

5.1.1. 只受自重的山体变形

由于山体在开挖前经过长期地形变化，地形变化后的山体如图4所示，需要考虑自重影响。

Figure 4. The von Mises stress in soil layer before excavation of the tunnel

图4. 隧道开挖前土层中的von Mises应力

5.1.2. 开挖单隧道山体变形

显示了隧道开挖后的应力分布，注意隧道周围有效应力的增加。

5.1.3. 开挖三隧道山体变形

在右半部山体中间再开挖一个隧道，进一步使山体发生复杂的应力和变形。

5.2. 模拟结果对比分析

数值模拟基于比较完善的工程力学理论，可以在一定程度上分析出并反映出隧道在施工过程中与地质环境相互作用的规律，对防止隧道开挖引发的次生灾害具有重要的参考价值。研究隧道开挖数量对土壤变形、应力分布的影响。

Figure 5. The von Mises stress in soil layer after excavation of single tunnel

图5. 单隧道开挖后土层中的von Mises应力

Figure 6. The von Mises stress in soil layer after excavation of three tunnels

图6. 三隧道开挖后土层中的von Mises应力

Figure 7. The plastic deformation in the zone near single tunnel after the excavation

图7. 开挖后单隧道附近区域的塑性变形

5.2.1. 只受重力作用下山体变化情况

这里需要单独考虑山体重力对自身的影响，是因为山体自身重力会引起土壤之间挤压产生内力，引起土体应力分布不均，不同部位发生沉降。由于假设条件中土体是均质的，所以处于同一水平截面位移变化情况相同。隧道开挖前由山顶到山底，土层中应力逐渐增大。

Figure 8. The surface subsidence of single tunnel

图8. 单隧道地面沉降图

5.2.2. 单隧道与多隧道开挖下山体变化情况

①应力分布情况(图5，图6)

隧道开挖前后所受应力处于不同的数量级，所以开挖隧道对土体所受应力分布有一定的影响。单隧道与多隧道开挖，都是在隧道与山体接触边界出所受应力较大、较集中，剩余部分所受应力较小。越靠近边界外部，应力越大；越靠近边界内部与山体接触，应力越小。

Figure 9. The horizontal displacement at the top surface of single tunnel

图9. 单隧道顶面水平位移图

Figure 10. The horizontal displacement at the top surface of three tunnels

图10. 三隧道顶面水平位移图

Figure 11. The plastic deformation in the zone near three tunnels after the excavation

图11. 三隧道开挖后隧道附近区域的塑性变形

Figure 12. The surface subsidence of three tunnels

图12. 三隧道地面沉降图

②隧道附近区域塑性变形分析(图7，图8)

单隧道与多隧道开挖引起的塑性变形相差不大，均是在隧道边界处发生的变形较大、较集中，剩余山体变形微小。变形范围都在0~1 mm。

③山体顶面水平位移变化(图9，图10)

单隧道开挖顶面水平位移变化比较规律，先下降后上升，最大变形约53 mm。多隧道开挖顶面水平位移变化出现两次波谷，先下降后上升，之后继续下降并高于之前下降的幅值，最后再上升。一次较大变形约为24 mm，最大变形约为62 mm。

④山体地面沉降变化(图11，图12)

开挖前是由于山体自身重力，使山体发生沉降。单隧道开挖由于隧道处于整个山体中心，所以隧道处发生的沉降较大约为104 mm，山体最两边发生沉降较小约为19 mm。多隧道开挖土体沉降发生不同变化，靠近中间隧道发生的沉降仍是最大的，约为176 mm；靠近两边隧道发生沉降也比其他地方明显，约为157 mm，小于单隧道开挖单独变形。

综上所述，隧道开挖所受应力分布及塑性变形与开挖隧道的位置相关，开挖隧道的数量与单个隧道开挖对土体的影响有相似性，但不是简单的线性叠加。若隧道开挖的体积比山体小的多，则隧道开挖的大小对山体变形影响较小。

6. 结语

本文应用COMSOL对山体隧道的开挖进行了数值模拟，采用弹塑性本构关系来描述山体的力学性质，并利用Drucker-Prager屈服准则，考虑隧洞开挖前后自重对山体受力达到稳态的影响。通过模拟山体开挖前后的拱顶和两侧山体沉降、变形，并将单隧道和多隧道开挖后山体的位移和应力变化情况进行对比，得出以下结论：

1) 隧道开挖过程中，由于考虑到自重及周围荷载对山体力学性质的影响，隧道上方及两侧位置的沉降和变形速率会随着时间的变化而变化，最后趋于稳定。

2) 隧道开挖应力主要集中在洞口周围，造成的土体变形主要位于隧道上方和左右两侧，变形主要以竖向地表沉降为主，关于隧道轴线对称两侧水平与竖直方向的变形规律相同。

3) 隧道开挖所受应力分布及塑性变形与开挖隧道的位置相关，多隧道开挖与单个隧道开挖对土体的影响有相似性，但不是简单的线性叠加，而是呈非线性的。隧道周围塑性区域的地表沉降和宽度是预测开挖过程中加固材料用量的重要参数，可为实际工程提供参考。