随着我国经济的飞速发展, 公路交通设施建设规模不断增长, 三车道以上的公路隧道设计逐渐成为主流^[1], 其中小净距隧道能够满足特殊地形、地质条件以及路线选择的限制, 逐渐得到广泛应用^[2]。

大断面小净距隧道受力特性与普通隧道不同, 较小的净距会导致双隧洞开挖对彼此的影响增大;较低的扁平率致使开挖后围岩和衬砌的力学特性发生明显变化, 拱脚处应力集中显著增大。因此, 无论从设计角度还是施工角度, 其结构受力状况都十分复杂^[3]。此外, 对于富水地区的隧道建设而言, 地下水也应成为设计施工过程中所要考虑的重要因素。孔隙水压作用会改变围岩的应力场和位移场, 围岩裂隙的发育又使得自身渗透性发生变化, 从而影响渗流运动, 流固耦合作用的不断发展严重影响了地下工程的安全稳定性^[4]。

以在建的大连市渤海大道一期工程大东山隧道为工程依托, 建立弹塑性应力-渗流-损伤耦合有限元计算模型;根据现场监测位移采用差异进化算法对损伤参数进行反演, 在所得参数的基础上, 对不同净距下, 隧道围岩的塑性区分布和不同水头高度下, 隧道围岩的损伤场演化规律进行分析;以隧道边缘观测点处的损伤值为指标, 研究了注浆圈渗透系数, 注浆圈厚度对损伤值和洞内涌水量的影响规律, 旨在为富水地区大断面小净距隧道工程的安全建设提供一定参考。

大连市大东山隧道全长2 223.3m, 全线为双洞双向八车道小净距隧道。工程所在地区夏季受太平洋副热带高压的控制, 多为东南季风, 气温较高, 降雨多, 降雨量约占全年降水量的73%。地勘资料表明, 拟建隧址区地下水埋深8.80～82.2m, 地下水位标高37.490～91.100m, 存在大量裂隙密集区域, 随着夏季降雨量的暴增以及原始地下水的存在, 洞内涌水量增大, 是工程安全的潜在隐患。

程序实现

本文引入考虑损伤的Drucker-Prager屈服准则对岩石塑性行为进行描述, 以有效应力原理和渗透系数演化方程为桥梁建立岩石弹塑性应力-渗流-损伤耦合模型, 并选取应力回映算法对非线性有限元问题进行求解。

计算过程中的基本假定: (1) 所研究的饱和岩体骨架为塑性损伤各向同性体, 并满足小变形假设; (2) 渗流为稳定流, 遵从Darcy定律; (3) 渗流液体为等体积等温的理想液体。

岩石微元体上的静力平衡方程为:

式中:f_i为体积力;σ_ij为总应力张量;Ω为问题求解域;n_i为边界法向余弦;t_j为作用在边界上的已知面力;Г₁为已知力边界;Г₂为已知位移边界;u_k^-为边界上的已知位移。

表征弹塑性体应变-应力关系本构方程为:

几何方程:

将式 (1) , (2) 和 (3) 联立, 即可求解位移模式下岩石弹塑性体内的应力应变状态。

将岩石损伤变量D表示为等效塑性应变ε^-p的函数^[5]:

式中:

式中:ε_p1, ε_p2, ε_p3分别为3个主塑性应变;等效塑性应变阈值ε^-p=0, 即等效塑性应变产生时有损伤演化;κ为试验所得正常数。

由Darcy定律可得不可压缩稳定渗流方程为:

式中:h为渗透水头;x, y, z为空间坐标;t为时间坐标;k_x, k_y, k_z分别为以x, y, z轴为主轴方向的渗透系数;S_s为单位贮存量;γ为水的重度;p为孔隙水压力。

水头边界和流量边界为:

式中:M₁为已知水头边界;f₁为已知水头边界值;M₂为已知流量边界;f₂为已知流量边界值。

在富水区岩体应力-渗流耦合场中, 一方面水体的渗流作用将改变围岩中的孔隙水压力, 从而改变围岩中的原始应力场状态;另一方面, 围岩应力-应变场的变化反过来影响围岩结构变化, 使得围岩的渗透性发生改变, 从而改变水体的渗流状态^[6]。用有效应力原理和岩石体积应变关系式能够对这一耦合过程进行描述。

由多孔介质的有效应力原理可知:

式中:σ'_ij为有效应力张量 (压为正、拉为负) ;p为孔隙水压力;α为等效孔隙水压系数, 0≤α≤1;δ_ij为Kroneker符号。

基于Kozeny-Carman公式的岩石渗透系数与体积应变间的关系表达式^[7]:

式中:n₀为初始孔隙率;ε_v为体积应变;K₀为初始渗透系数。

塑性损伤耦合主要通过势函数、加载函数和一致性条件进行相互影响。采用Druker-Prager屈服准则对岩石塑性应力、变形特征进行描述:

式中:p为静水压力;J₂=s_ijs_ij, s_ij为偏应力;s=σ-p (σ) I;, i, j=1, 2, 3;为考虑损伤作用的黏聚力;D为损伤变量。

式中:为初始黏聚力;c_r为岩石破坏时的黏聚力;ζ为材料损伤参数, 取值范围为0≤ζ≤1。

本文采用固体和渗流有限元理论对应力场平衡方程 (1) ～ (3) 和渗流场平衡方程 (6) 进行离散得:

式中:[K_m]为应力场总刚度矩阵;{u}为节点位移列向量;Δu为增量位移;{f_m}为节点力向量, 包括体力、面力和孔压的等效荷载;[K_s]为渗透矩阵;{h}为水头列向量;{f_s}为自由项列向量。

本文利用牛顿-拉弗森法对式 (12) 进行迭代求解得出位移增量, 通过位移增量计算得到应变增量, 进而得到应力增量。为避免非法应力的存在 (计算所得应力脱离屈服面) , 本构积分算法采用隐式返回映射算法, 对计算所得应力进行塑性修正 (见图1) 。

返回到屈服面的光滑圆锥面时的一致切线模量为:

式中:I_d是体积张量, a, b, c, d是系数。

返回尖点处时的一致切线模量为:

进行耦合计算时, 先利用弹塑性损伤程序计算出岩石体积应变, 再按照式 (9) 计算各单元的渗透系数矩阵并带入方程 (12b) 中进行渗流场求解, 将得到的孔隙水压通过有效应力原理式 (8) 对应力状态进行修正, 反复迭代直至前后两次求得的应力场、渗流场都满足收敛条件为止。具体有限元程序的实现见文献[8]。

本文选取了全线净距最小的断面K06+940 (净距17m, 单洞外轮廓跨度16m) 作为研究对象, 建立隧道平面应变有限元计算模型。模型高为50m, 宽为100m。采用四边形网格对围岩体进行剖分, 共划分777个单元和859个节点。模型两侧施加水平方向的位移约束, 底面施加固定位移约束。参照地质勘察报告, 此断面处岩体质量等级为Ⅳ级, 其基本计算参数选取如表1所示。有限元计算模型网格划分如图2所示。

由于围岩结构体的复杂性和非确定性, 参数选取问题很大程度上影响了数值计算结果的稳定性。因为隧道的开挖力学行为是围岩体力学性质的真实反映, 所以利用监测信息能够获得合理的围岩体力学参数, 这种方法即为反分析法, 也称参数识别。反分析法的本质是优化问题。差异进化算法是一种新型全局优化算法, 对初始值无要求, 收敛速度快, 能够适应各种非线性函数, 具有传统遗传算法所没有的优势^[9]。

本文结合建立的有限元正算模型, 编制基于差异进化算法 (DE) 的损伤参数耦合反演程序, 程序计算流程如图3所示^[10]。

现场位移监测点布置如图4所示。根据现场监测数, AB=10.388mm, AC=9.370mm, AD=8.254mm, AE=8.937mm, BC=0.566mm, DE=0.788mm初步给定围岩损伤参数, 明显损伤时黏聚力c_r=30kPa, 损伤参数ζ=0.7, κ=1 000, 代入图3所示程序进行参数反演计算, 所得结果如表2所示。

隧道开挖会引起围岩应力重分布, 应力作用下围岩的弹塑性区域就会发生变化。当双洞距离足够远时, 开挖引起的应力场互不干扰;随着双洞空间距离的缩小, 相互影响范围发生相交或重叠, 使得中间岩柱区域出现塑性区贯通, 进入破坏状态, 可能引起隧道整体失稳, 从而影响围岩与结构的安全问题。

在不改变其他参数, 不引入渗流方程的条件下, 对不同净距下 (17, 15, 13m) 围岩塑性区的分布进行数值模拟计算, 计算结果如图5所示。

由图5可知:隧道净距对围岩塑性区分布的影响是显著的, 随着净距的减小, 隧道中间岩柱的塑性区逐渐增大, 隧道外侧塑性区增大相对较小。在净距为13m时, 围岩塑性区在中间岩柱部位已经贯通, 围岩应力已经呈现较差状态, 严重影响到施工安全问题。数值结果也表明, 从围岩稳定性角度而言, 在当前围岩等级下设计所给的17m净距较为合理。

在当前设计净距下, 加设渗流方程, 对流固耦合作用下围岩损伤区进行数值模拟研究。以计算模型下表面为基准, 分别设置水头高度H为15, 20, 25m, 在模型左右两侧施加沿重力方向呈梯度变化的水头压力。取围岩渗透系数K_x=K_y=7.8×10^-5cm/s, 初始孔隙度e₀=2.7×10^-3。不同水头高度下孔隙水压力如图6所示。

由图6可知, 隧道开挖后, 孔隙水由高势能面流向低势能面, 不断向隧道内部进行渗透。在不同水头高度下, 围岩损伤区分布情况如图7所示。

计算结果表明, 围岩损伤区主要分布在隧道两腰处并随着地下水水头高度的增加不断向两侧扩张, 甚至有在中间岩柱发生贯通的趋势, 地下水的渗流作用加剧了小净距隧道开挖后围岩损伤区的演化, 对工程安全造成了严重隐患。

在地下水水头高度为45m的情况 (见图7c) , 对周围围岩进行注浆, 注浆参数如表3所示。分别在计算模型左边开挖洞的拱顶、拱肩和拱底设置损伤值观测点如图8所示。

设置注浆圈厚度为4m并保持不变, 成倍减小注浆圈系数, 记录不同注浆圈渗透系数下各监测点损伤值变化规律如图9所示。

由图9可知, 在注浆圈厚度保持不变的情况下, 随着注浆圈渗透系数降低, 隧道观测点处的损伤值均有不同程度降低, 注浆圈的存在抑制了围岩中地下水的渗流作用, 减小了孔隙水压对围岩应力场的作用。当围岩渗透系数与注浆圈渗透系数比n>40时, 损伤值变化趋势减小。

图10给出了在注浆圈渗透系数不变的情况下, 增大注浆圈厚度时, 拱肩、拱脚和拱腰1号点位置处损伤变化规律。

由图10可知, 隧道边缘观测点处的损伤值随着注浆圈厚度的增加不断减小, 当注浆圈厚度t>8m后, 损伤值变化趋势减小, 接近为0。

对不同注浆圈厚度下, 单个洞室内的涌水量Q进行计算, 结果如图11所示。

由图11可知, 涌水量Q随着注浆圈厚度的增加逐渐减小, 当注浆圈厚度从2m增大到6m时, 涌水量Q由11.95m³/ (m·d) 减小到5.01 m³/ (m·d) , 减小趋势十分明显;当注浆圈厚度>8m时, 涌水量变化趋势几乎为0, 增加注浆圈厚度已经无法有效控制涌水量大小。

由上述分析可知, 在本文计算条件下, 当n>40时且t>8m时, 减小注浆圈或增大注浆圈厚度对减小隧洞周边孔隙水压, 抑制隧洞围岩损伤演化的作用已不明显。因此, 注浆圈渗透系数和注浆圈厚度的选取存在经济上的最优解。

1) 以在建大东山隧道为工程依托, 建立弹塑性应力-渗流-损伤耦合有限元计算模型, 利用差异进化算法, 通过现场位移监测值求解损伤参数。

2) 围岩塑性区随着隧道净距的缩小不断扩张, 并在中间岩柱部分有贯通现象。数值计算结果表明:在当前围岩等级下, 隧道净距设计较为合理。

3) 围岩损伤场随着地下水头高度的增加演化加剧, 地下水对围岩稳定性存在较大影响。

4) 注浆圈的存在能够有效抑制地下水渗流作用, 减少隧洞边缘观测点的损伤值。