隧道力学特征及数值模拟方法

2 隧道力学特征及数值模拟方法
2.1 隧道开挖生成的围岩二次应力场特征
岩体在开挖前处于初始应力状态，初始应力主要是由于岩体的自重和地质构造所引起的。

在岩体进行开挖后改变了岩体的初始应力状态，使岩体中的应力状态重新分布，引起岩体变形甚至破坏。

在这个时间工程中，地层应力是连续变化的，特别地，洞室开挖后在未加支护的情况下，地层应力所达到的新的相对平衡称为围岩的二次应力状态。

一般来说，二次应力场是三维场。

在隧道施工过程中，横向的二次应力作用使得洞周围岩的应力状态和变形状态发生了显著的变化，可将洞周围岩从周边开始逐渐向深部分为4 个区域：（1）松动区由于施工扰动（例如施工爆破），区内岩体被裂隙切割，越靠近洞室周围越严重，其内聚力趋近于零，内摩擦角也有所降低，强度明显削弱，基本无承载能力，在重力的作用下，产生作用在支护上的松动压力。

（2）塑性强化区这一区域是围岩产生变形的根源。

隧道开挖后破坏了地层的原状力线，在洞体四周产生了很高的应力集中，此时，该处只存在切向应力和指向隧道中心的径向不平衡力，切向应力由承载拱承担，而对于径向应力，毛洞是无法承担的，所以要释放（在有支护的情况下一部分被初期支护承担）。

这就造成了洞体开挖后四周的围岩向隧道中心发生位移，周边的径向应力逐渐趋向零，而切向应力随着径向位移而增大。

这一应力状态的变化导致岩体从初始的二轴（这里只考察平面应力状态）受压状态转变为单轴受压状态，使得这一区域围岩处于非常不利的受力状态，当这一应力状态超过岩体的强度极限时，洞室周围出现了塑性区域或者破坏区域，产生塑性变形。

如果洞室周围塑性区域扩展不大，随着径向位移的出现，地层塑性区域达到稳定的平衡状态，围岩没有达到承载能力的极限值；但是如果塑性区域继续扩展，则必须采取支护措施约束地层运动，才能保持洞室围岩处于稳定状态，这时为了阻止地层运动，就显出塑性变形压力。

(3) 弹性变形区域这一区域内岩体在二次应力作用下仍处于弹性变形状态，各点的应力都超过原岩的应力，应力解除后能恢复到原岩应力状态。

其次，开挖面前方地层对已开挖区域的围岩有某种程度上的纵向支撑作用，即产生纵向的承载拱，承载拱的跨度约为一倍洞径。

所以
在开挖过程中不同横断面上的二次应力分布和变形是不同的，这种开挖面效应使得在开挖面前方和后方一定范围内的二次应力场呈三维分布状态，这种空间效应可以通过隧道开挖后的位移场反映出来(如图2.1)。

图2. 1隧道开挖后的位移场
(4) 原岩状态区岩体未受影响，仍处于原岩状态。

新奥法认为围岩压力类型可以相互转化，即变形压力可以转化为松动压力。

支护结构对围岩起到约束作用，如果支护强度不足，围岩位移过大进入松弛变形阶段, 则形成了松弛压力。

松弛压力的实质是部分围岩脱离母体后作用于支护结构上的重力。

围岩中常分布着大量的裂隙和结构面，而这些薄弱面的剪切强度远低于岩石块体强度，而且不能承受拉力的作用，因此岩体力学性质往往表现为各向异性，在特定方向上的强度远低于岩石的块体强度。

当结构面和隧道临空面形成不利组合时，失去支撑的围岩块体脱离母体(图2.2)，其重力作用于支护结构上形成松弛压力，这时围岩压力将不在取决于围岩的特征曲线。

图2.2围岩块体塌落是松动压力的主要来源
2.2隧道施工过程中围岩对支护作用力的探讨
隧道施工时，由于承载拱效应，原始的地层应力并不是全部转化为作用在结构上的荷载，即使在隧道建成若干年后，作用在隧道衬砌上的压力仍然小于初始应力。

这已被大量的工程测试资料所验证。

对于这种现象，比较权威的解释是[15] [15]王梦恕•地下工程浅埋暗挖技术通论[M] •安徽：安徽教育出版社，2005，隧道开挖后洞室周围地层应力的释放，隧道的拱形形状以及地层内部摩擦力等导致承载拱发挥作用，周围地层应力进行重分布产生两种变化，即一部分被释放，另一部分想深部和其他方向转移。

当施作衬砌支护后，地层应力的释放过程受到拟制，一部分释放荷载作用在衬砌结构上，这部分荷载的大小正是我们需要了解的作用于衬砌结构上的压力。

围岩的松动塌方与提供支护的时机有关，如果支护愈早，提供的抗力就愈大（图2.3），围岩就能稳定。

反之，支护迟，提供的支护抗力愈小，不足以维持围岩的稳定，松动区中的岩体在重力作用下会松动塌落，所以要维持围岩稳定，既要维持围岩的极限平衡，还要维持松动区内滑移体的重力平衡。

如果为维持滑移体重力平衡
o -------- Advancing tunnel
图2.3与距掌子面位置相关的支护压力[16] [16]刘波，韩彦辉.FLAC原理、实例与
应用指南[M] •北京：人民交通出版社，2005
所需的支护抗力小于维持围岩极限平衡状态所需的支护抗力，那么只需要松动区还保持在极限平衡状态之中，松动区内滑移体就不会松动塌落。

反之，则会松动塌落。

由此，我们可以把维持松动区内滑移体平衡所需的抗力等于维持极限平衡状态的抗力，作为围岩出现松动塌落和确定P imin条件。

要确定最佳支护结构或最佳支护时间，必须确定最小围岩压力P imin。

最小围岩压力和围岩允许最大位移两者是等价的（图2.4）。

目前，无论是确定最小围岩压力还是确定围岩允许最大位移都没有好的计算方法。

对于一1的情况，可以用下面方法估算[14] [14]何满潮，黄润秋王金安等.工程地质数值方法[M] •北京：科学出版杜，2006:
图2.4围岩特征曲线和支护特征曲线［16］［佝刘波，韩彦辉.FLAC 原理、实例与应
用指南［M ］ •北京：人民交通出版社，2005
由岩体力学可知，在• =1的情况下，围岩松动区内的滑裂面为一对对数螺线。

假设 松动区内强度已大大下降，可以认为滑移岩体己无自承作用，以致于松动区内滑移 体的全部重量都要由支护抗力 P imin —来承受，由此有：
P i min =G (2.1)
考虑到实际情况，真正作用在支护结构上的压力应当是重力与变形压力的叠加，
则式(2.1)应该写为：
P i min 二 2 G 式中：G 为滑移体的重量；b 为滑移体的地宽。

滑移体的重量可以近似取下式：
'b ( R max - r 0 ) G 二 ------------
2 AneesMJdEddns
(2.2) (2.3)
frrwvrd radial diBpfaaQorrwrrt u*
式中：R max为与P imin相应的允许最大的松动取半径；为岩体容重
(2.4)
按塑性区半径公式即可得:
计算R max 时，采用的C 值应该再降低
P imi n 的大小主要取决于松动区半径R max ，当原岩应力愈大，C ,「值愈低和C ，值 损失愈多时，则R max 和P imi n 就愈大。

此外，还与岩体构造状况，施工爆破情况、外 界条件等有关，这些都会影响围岩 c 值的降低。

由于隧道围岩力学性质的复杂性及不确定性，要纯粹利用解析的方法准确计算 围岩对支护的压力是很困难的，但是可以通过数值模拟及实测数据预测围岩对支护 的作用。

来三维力学分析结果及国内外大量的实测数据均表明，原始地层应力的释 放率与地面沉降和拱顶下沉之比有很好的一致性 [14] [14]何满潮，黄润秋王金安 等•工程地质数值方法[M] •北京：科学出版杜，2006,即 = U 地面/U 拱顶 (2.6)
而 / ~0 (2.7)
式中匚v 为实际作用在结构上的压力，匚0为原始地层应力
于是
:二v = ：•匚 =（U 地面/U 拱顶）- (2.8) 对于拱部：
:-v - (U 地面/ U 拱顶)H
(2.9)
对于边墙：
j v — (U 地面/ U 拱顶)K 0 h
(2.10)
对于底部：
■
w
v - (U 地面/ U 拱顶)H (2.11)
DL 上式中，w 为衬砌结构的自重；H 为覆盖土厚度；D 为隧道结构跨度；i 为隧道侧 向任意一点力地面的距离，L 为闭合长度。

对于作用在结构上的围岩压力与覆跨比的关系，由上面各式可以看出，作用在 结构上的围岩压力与荷载释放率 6密切相关，而艿即为地面沉降与拱顶下沉之比， 因此，只要研究U 地面/U 拱顶与
-sin 半 YI I 1 +sin J 1 _sin :
2 sin : j (2.5)
R max
P +cctg ④ i min +浊®丿
H / D的关系，就可以弄清楚作用在结构上的围岩压力与覆跨比的关系了。

日本的岛田隆夫通过试验得到以下公式：
U地面/U拱顶》e」H/D)(2.12) 式中，H / D为覆跨比，〉、■:是与施工方法有关的常数。

Atkinson等人根据模型试验得到下式；
U 地面/U 拱顶=1.0 - ： H / D (2.13) 将2.12式和2.13式代入2.6式和2.11式就可以算得隧道在开挖过程中的围岩释放率和围岩作用在支护上的作用力。

2.3数值模拟方法
随着计算机技术的迅速发展，借助数值模拟方法与计算机图形、图象技术、可视化技术相结合，对地下工程的开挖步骤、支护工艺的工程性态(稳定和变形)模拟
和过程再现已成为现实，这种工程性态模拟和过程的图形显示技术被称为计算机仿真技术。

在数值模拟中，不仅仅可以模拟岩体中的断层、节理、裂隙等地质结构面，而且还可以模拟分布开挖、支护等施工过程，揭示不同施工步骤、施工工艺以及不同支护条件下的应力与位移。

尤其采用不连续分析技术，模拟和显示洞室开挖过程，洞顶及边墙有些部位岩块失稳而下落或滑移，为支护设计提过可靠依据。

这些都是解析方法难以实现的。

本文利用岩土工程界最为常用的有限单元法和有限差分法，对厦门海底隧道不同施工段进行了施工力学行为模拟分析。

2.3.1弹塑性有限元分析
岩体弹塑性本构关系是岩石主要非线性问题之一。

岩石的弹塑性是指岩石材料的应力应变关系在屈服之前呈线性关系，当应力达到屈服应力时，应力应变关系就变为非线性。

由于弹塑性模型中应变不仅依赖于受载的应力状态，而且与加载路径有关，因此一般弹塑性本构关系不能用应力应变全量关系准确描述，只能用能反映与加载路径有关的应力应变增量关系描述。

在岩石非线性本构关系有限元分析中，一般采用初应力法和初应变法求解非线性平衡方程组。

初应力法是将荷载以微小增量形式逐级加在模型上，每加一级荷载增量：dF i /，就会产生相
应的位移增量d j ?＞应变增量〔d ;i［和应力增量J o对于具有初应力的弹塑性应力应变本构关系可以写成：
- b p td ：(2.14) :d；「o ： - - b p id (2.15) 其中D p 1为塑性矩阵，它与加载前的应力水平有关，而与应力增量无关。

初应力法是通过对；0打的处理将应力修正到正确的水平上，初应力心亠？不仅与加载增量前应力水平有关，还与本级所加荷载增量引起的应变增量％；［有关。

增
量形式
平衡方程为：
K 0 Id二加.；：，dF 1 (2.16)
式2.16中，K。

为线弹性计算中总刚度矩阵；】为校正荷载项，由式2.17 决定：
：dF 1 - 7 I B：b Ip Id ；：dv (2.17) 由于0F ?随位移变化而变化，所以计算时必须进行迭代求解。

初应力法求解按照以下步骤实现：
(1) 把全部荷载划分成若干个增量，在每一级增量段内，按照增量弹塑性平衡方程进行求
解。

(2) 计算各单元的应力增量及当前应力
仏}j = b Id 3 }j
丿{d j }j = b ］{d® }j (2.18)
匕片=匕L +{d§片
下标i表示第i级荷载增量；j表示第j次迭代。

(3) 根据岩石的屈服准则，由各单元应力判断单元是否屈服，对于塑性单元，计
算应力修正项并修正应力
(2.19)
(4) 塑性单元通过修正项^cp ［计算等效节点力，所有塑性单元的等效结点力叠
加构成总的修正荷载矢量
I dFi ] =X t d CT p }. dV (5) 在修正荷载作用下进行下次迭代运算，此时基本方程为
［K ］d.j =「dF i 重复进行(2)〜(5)步计算直至所有的塑性单元达到收敛精度要求。

在进行下一步 的荷载增量计算
(6) 重新施加下一级荷载增量fdF i.J ，重复计算
(1)〜(5)步，直至计算完毕。

通过 累加各级载荷作用的计算结果，求得总位移 和总应力匚门。

一般初应力法的收敛速度比较缓慢，因此通常采用常刚度和变刚度法相结合的 方法加速收敛。

在ANSYS 中岩土工程问题的分析一般使用 Drucker--Prager 屈服模型(图2.5)， 其等效应力表达式为：
(2.22)
其中：二m - ；「x ，「z /3为平均应力或静水压；'S ；■为偏应力；1为材料常
数；［M ］为mises 屈服准则中的［M ］(2.20) (2.21)
图2.5德鲁克一普拉格屈服准则
上面的屈服准则是一种修正的 Mises 屈服准则，他考虑了静水应力分量的影响,
静水应力（侧向压力）越高，则屈服强度越大，材料常数 1的表达式如下：
p
2 sin $
*3（3 _ sin © ）
屈服准则的表达式如下：
6C cos
-・y
<3（3 -sin © ）
最后的屈服准则表达式为:
(2.25)
2・3・2有限差分法分析
⑴概述
有限差分法是较早用于求解给定初值和（或）边值微分方程组的数值方法。

在有 限差分法中，基本方程组和边界条件（一般为微分方程）近似地改用差分方程（代数方 程）来表示，及由空间离散点处的场变量（应力、位移）的代数表达式代替。

这些变量 在单元内是非确定的，从而把求解微分方程的问题转化成求解代数方程的问题。

而有限元是将连续求解域离散成一组有限个相互连接的单元体，利用每一个单
-a
(2.23)
(2.24)
元内假设的近似函数来分片求解区域上待求的未知场函数。

单元内的近似函数通常
由未知场函数及其导数在各节点的数值和插值函数来表示，这样有限元分析中的未
知函数及其导数在各个节点上的数值就成为新的未知量，一经求出这些未知量，就可以通过插值函数计算出各单元内场函数的近似值，进而得到整个求解域上的近似解。

有限差分法和有限元法都产生一组待解方程组。

尽管这些方程是通过不同方式推导出来的，但是两者产生的方程是一致的。

另外，有限元程序通常要将单元矩阵组合成大型整体刚度矩阵，而有限差分则无需如此，因为它相对高效的在每个计算步重新生成有限差分方程。

在有限元法中，常采用隐式、矩阵求解方法，而有限差分法则通常采用显式、时间递步法解代数方程。

FLAC 采用拉格朗日分析方法，由于它不需要形成整体刚度矩阵，大变形计算
时在每个计算步都很容易修正坐标。

位移增量施加到坐标上以致网格随着材料发生移动和变形，这就是所谓的拉格朗日法；若材料移动和变形是相对于固定的网格，就是与此对应的欧拉法。

(2) 有限差分法基本方程
在弹性体上用相隔等间距h 并平行于坐标轴的两组平行线划分成网格(图2.6)
其代入式2.26,得
假定h 充分小，因而可以不计它的三次幕及更高次幕的各项，则式 可以简化为:
设f =f x, y 为弹性体内某个连续函数，它可能是某个应力分量或位移分量，也可 能就是应力函数、温度、渗流等。

这些函数，在平行于某轴的一根格线上, 坐标的变化而改变。

在临近节点 0处，函数f 可以展开为泰勒级数：
只随
…。

to
1
X
-X o 匚
j 2 f -2
::X
’ 2
1 X _Xo ) +_
3! 丿0
在节点3和节点1，x 分别等于X o
(X _X 。

j +…
o
-h 和X o - h ，即x o -h 分别等于
(2.26)
-h 和h ，将
f 3 二 f o -h
f l = fo
h
h 2
+ —
2
J f
;x 2 g 2
f
r 2
g 丿o
(2.27)
(2.28)
(2.27)(2.28)
联立求解(2.29)(2.30),可得到差分公式:
f
1 ■ f3 ~-
2 f0
同样，可以得到
:f _ f2 - f4
：：x 0 2h
f2 J - 2 f°
h2(2.29) (2.30) (2.31) (2.32) (2.33) (2.34)
2
r 2
2丿0
公式(2.31)(2.34)是基本差分公式，通过这些公式可以推导出其他的差分公式
例如，利用公式(2.31)(2.33)，可以导出混合二阶导数的差分公式
(2.35)
用同样的方法，由公式(2.32)(2.33)可以导出四阶导数的差分公式。

(3) 平面的有限差分拉格朗日法数值原理
对于平面问题，将具体的计算对象划分成四边形网格域的有限差分网格，每个
图2.7有限差分单元划分示意图
三角形网格域的有限差分公式可以用高斯散度定理的广义形式推导得出:
在面积A 上定义f 的梯度平均值为:
边的平均值。

平面问题有限差分法是基于物体运动与平衡的基本规律。

最简单的例子是物体
其中：
s
为绕闭合面积边界积分，n j 为对应表面s 的单位法向量，f 为标量、
s “i
fds 二 A
工dA
矢量或张量，X j 为位置矢量，ds 为微量弧长，A 表示对整个面积A 积分
势⑴
(3.36)
■:f -：
A A 丄 dA A ;X i
将式(2.36)代入上式可得:
对于一个三角形的网格域，式(2.38)的有限差分形式为:
其中，As 为三角形网格域的边长，求和是对该三角形的三个边进行;
(3.37)
(3.38)
(3.39) ■； ＜取
该
1 A -n j dS
二 一 ' f n i As
A s
质量(m)、加速度(du/dt )与施加力F 随时间变化的关系，牛顿定律描述的运动方程 为：
du m ——=F dt
图2.8物体对随时间变化作用力的响应
当几个力同时作用于该物体时，如果加速度趋于零，即
a F = 0 (对所有作用
力求和)。

式(2.40)也表示该系统处于静力平衡状态。

对于连续固体介质，式 (2.40河
以写成如下广义形式：
dt
式中，T 为物体的质量密度，t 为时间，薯为坐标矢量分量, 力)分量，5为应力张量分量 禾I 」用式(2.39)，将f 替换成网格域每边平均速度矢量，这样，网格域的应变速率 ej 可以用网格点速度的形式表述：
式中，(a)(b)是三角行边界上两个连续的网格点
如果节点间的速度按线性变化，式(2.42)平均值与精确积分是一致的。

通过式
(2.42)、(2.43河以求出应变张量的所有分量
(2.40)
► F(t)
(2.41)
晶为重力加速度(体
e
ij
(2.42)
.:u i
j
-(b)
2A S
n j As
(2.43)
P
du
:x j
1
;:u j
1 '(a ) 2
根据力学本构定律，可以由应变速率张量获得新的应力张量：
G =M （G ，e ij ,k ）
（2.4
4）
式中，M … 表示本构定律的函数形式；k 为历史参数，取决于特殊本构关系； 表示“由……替换”。

在FLAC 程序中，单元应变率是计算各主要参数的纽带，由于非线性应力应变
本构关系不具有唯一性，一般用增量形式表示，当已知单元旧的应力张量和应变速
率（应变增量）时，可以通过式（2.44）确定新的应力张量。

例如各向同性最简单的弹性 本构关系张量差分形式如下：
式中At 为时间步，G ， K 分别是剪切模量和体积模量。

在一个时步内，单元的 有限转动对单元应力张量有一定的影响。

对于固定参照系，此转动使应力分量有如
G ：二 5
丸
在大变形计算中，先通过式（2.46）进行应力校正，然后利用式（2.45）或（2.44）计算
等前时步的应力。

一旦计算出作用在网格域上的应力后，就可以确定作用到每个网 格点上的等价力。

在每个三角形子网格域中的应力如同在三角形边上的作用力，每 个作用力等价于作用在相应边端点上的两个相等的力。

每个角点受到两个力的作用, 分别来自各相邻的边。

因此:
F j =丄6（n （ S （1）+ n f S （2））
（2.48
）
由于每个四边形网格域有两组两个三角形域，在每组中对每个角点处相遇的三
角形节点力求和，然后将来在这两组的力进行平均，得到作用在该四边形网格点上 的的力。

在每个网格点处，对于所有围绕该网格点四边形的力求和
x
F j ,得到作
用于该节点的纯粹节点力矢量。

该矢量包括所有施加的荷载作用以及重力引起的体 力F g
i ，有：
Gj
：二
CT j + < 6j
ij
2 |. -
K …—G e kk • 2 G e j • = t
3
(2.45)
(2.46)
其中，■ ■打 =-
；：U i ：：u j
(2.47)
' F j ^g i m g (2.49) 其中m g是凝聚在网格点处的质量，定义为连接该网格点的所有三角形网格域质量和的三分之一。

如果四边形网格域不存在，例如空单元，则忽略对F i的作用；如果物体处于平衡状态，或处于稳定的流动状态，在该网格点处的-F i被视为零。

否则根据牛顿第二定律的有限差分形式，该网格i将被加速：
5宀、F j* (2.50)
m
其中，上标表示确定相应变量的时刻。

对于大变形问题，将式(2.49)再次积分，可以确定出新的网格点坐标：
x「弋叮⑴耳(2.51)
注意到式(2.50)和(2.51)都是在时间段中间，所有对中间差分公式的一阶误差项消失，速度产生的时刻，与节点位移和节点力在时间上错开半个时步。

(4) 三维有限差分拉格朗日法数值原理
对于三维问题，先将具体的计算对象划分成六面体有限差分网格，每个离散化后的立方体网格可进一步划分出5个常应变三角棱锥体子网格。

每个四面体网格点
是从1到4标记，约定网格点号与所对的面同号
node*
1
图2. 9立方体网格划分成5个常应变三角棱锥体网格应用高斯散度定理于三角棱锥形网格，可写成:
V v
i,j dv = j s v i n j dS
(2.52)
式中积分分别是对棱锥体体积和面积进行的，［n］是锥体表面的外法线矢量。

对于恒应变速率棱锥体，速度场是线性的，并且［n］在同一表面上是常数。

因此，通过对式(2.52)积分，可得：
4
VV i,j =送v(f n j(f S C ) (2.53)
f J
式中的［f］表示与表面［f］上的附变量相对应，V是速度分量f的平均值。

对于线性速率变分，有
_ (f) 1 I
V i 二一' V i (2.54)
3
I zfcl 尹
式中的上标，表不关于网格点，的值。

将式(2.54)代入(2.53河得到网格点和整个网格域的关系：
1 4 | °彳
VV j,j = - £v i O J n j(S(f) (2.55)
3 i _i i ±i 二f
如果将式(2.52)中的V i用1替换，应用发散定律，可得：
4
n j f S f =0
f 土
利用上式并用V除以式(2.55河得:
应变速率张量的分量可以表述成:
由此可进一步推导网格点的运动方程，同样基于物体运动与平衡的基本方程。

(2.56)
v i,j
1
3V
(2.57)
ij 6V
4
送(v：n j O+v j n i()S©)
l T
(2.58)。