CN114676486A - 江水渗流对穿江隧道开挖影响的分析方法 - Google Patents

Abstract

本发明公开了江水渗流对穿江隧道开挖影响的分析方法，建立隧道开挖力学分析模型和开挖隧道江水渗流分析模型；基于隧道开挖力学分析模型，进行隧道开挖过程的力学分析，得到隧道开挖面上的开挖荷载；基于开挖隧道江水渗流分析模型，进行渗流场分析和开挖荷载作用下的隧道应力与变形分析，得到洞周围岩内部微结构的变化；根据岩体内部微结构变化，采用岩体渗透性演化方程进行洞周围岩渗透张量的计算，并将其传递到渗流分析模型中；采用双场交叉迭代算法进行江水渗流过程开挖隧道围岩应力、变形与渗流过程的求解。本发明还分析江水变化与岩体结构特征变化在穿江隧道开挖施工中对江底沉降变形的影响，具有创新意义，工程应用前景广泛。

Description

江水渗流对穿江隧道开挖影响的分析方法

技术领域

本发明涉及数值模拟分析技术领域，尤其是江水渗流对穿江隧道开挖影响的分析方法。

背景技术

进入21世纪以来，随着国民经济的高速发展，城市化进程不断加快，人口越来越向大城市集中，从而要求城市必须不断扩大其规模以适应城市人口的急剧增加。因此，为满足人们对出行、生活等方面越来越高的需求，我国对地下工程和隧道工程发展的步伐越来越快，而其中穿江隧道与桥梁相比，具有不受气象条件制约、不影响航运、防灾减灾性能好、战略意义高等诸多优点。

对于大型过江隧道，其在修建过程面临开挖隧道直径大、含水量高、自稳性差等重要问题，施工过程中对周边围岩扰动大，易出现地表变形和位移过大等情况，可能导致开挖面塌陷、江水突涌等安全事故。富水地层隧道修建，不可避免的遇到与水有关的岩体水力耦合问题。受岩体结构及其所处工程地质环境的复杂与不确定性影响，水力耦合问题研究一直极具挑战。

针对隧道开挖扰动响应的模拟，目前研究往往主要考虑力学过程，较少考虑岩体结构特征与渗透特性演化等的影响，使得研究有待进一步深入。江水渗流与岩体结构作为影响隧道开挖变形的重要因素，在实际穿江隧道设计与施工过程得到广泛关注。因此，在穿江隧道开挖扰动中考虑江水渗流影响效应的研究，不仅对深化涉水隧道水力耦合研究具有参考意义，同时还具有十分重要的工程意义。

发明内容

为了克服上述现有技术中的缺陷，本发明提供江水渗流对穿江隧道开挖影响的分析方法，基于岩体水力耦合分析方法，对江水渗流过程与围岩损伤变形过程进行耦合，对穿江隧道开挖过程中江水渗流的影响效应进行了***的分析，应用前景良好。

为实现上述目的，本发明采用以下技术方案，包括：

江水渗流对穿江隧道开挖影响的分析方法，包括以下步骤：

S1，建立隧道开挖力学分析模型；

S2，基于隧道开挖力学分析模型，通过将隧道开挖区域挖除，建立开挖隧道江水渗流分析模型；

S3，基于隧道开挖力学分析模型，进行隧道开挖过程的力学分析，得到隧道开挖面上的开挖荷载；

S4，基于开挖隧道江水渗流分析模型，进行初始渗流场分析，分析开挖隧道江水渗流分析模型中的水压力；将开挖隧道江水渗流分析模型中水压力的分析结果代入隧道开挖力学分析模型中，基于隧道开挖力学分析模型，进行开挖载荷和水压力共同作用下的隧道应力与变形分析，分析洞周围岩内部微结构的变化；

S5，建立岩体渗透率演化方程，并根据洞周围岩内部微结构的变化，计算洞周围岩的岩体渗透率；

S6，将岩体渗透率代入开挖隧道江水渗流分析模型中，继续进行渗流场分析，更新开挖隧道江水渗流分析模型中的水压力，按照步骤S4～S5的方式，采用交叉迭代算法，分析江水渗流过程中隧道的应力与变形，以及洞周围岩内部微结构的变化。

步骤S1中，建立隧道开挖力学分析模型，具体包括以下步骤：

S11，根据穿江隧道工程地质条件、隧道尺寸结构与布置以及江水条件，采用八节点六面体单元，建立隧道初始有限元计算模型即初始模型；其中，初始模型的顶部为江底，整个模型均处于江水位以下，Z方向为竖直方向，围岩为饱和状态；

S12，对初始模型的底部和两侧均施加零位移法向约束；

S13，根据实际隧道工程特点，采用初始自重应力计算或采用侧压力系数分析的方式，设置初始模型的初始应力；

所述初始自重应力计算的具体方式为：根据江水深度在初始模型顶部设置均匀分布的水压力荷载，在隧道自重荷载作用下对初始模型的应力进行计算，根据应力计算结果设置初始模型的初始应力；

所述侧压力系数分析的具体方式为：根据侧压力系数K_c，采用如下计算公式，对初始模型的的初始应力进行设置：

σ_z0＝γ_wH_w+γ_rH_r

σ_x0＝σ_y0＝K_cσ_z0

其中，σ_z0为竖直方向的初始应力，σ_x0为水平面沿x方向的初始应力，σ_y0为水平面沿y方向的初始应力，γ_w和γ_r分别为水和围岩的重度，H_w和H_r分别为江水深度和围岩埋深。

步骤S2中，建立开挖隧道江水渗流分析模型，具体包括以下步骤：

S21，在隧道开挖力学分析模型的基础上，将隧道开挖区域挖除，建立隧道江水渗流有限元计算模型即江水渗流模型；

S22，设置江水渗流模型的边界与初始条件；

江水渗流模型的边界具体为：模型两侧和底部边界设为隔水边界；顶部设为固定水压力边界，顶部的固定水压力边界值根据江水位高程确定；在隧道开挖后，隧道洞周边界水压力降为1标准大气压；

江水渗流模型的初始条件具体为：初始时，江水渗流模型的内部水压力从顶部到底部呈梯形分布。

步骤S3中，通过在隧道开挖面节点上施加等效的开挖荷载来模拟隧道施工开挖过程；

隧道开挖面节点上所施加的开挖荷载的有限元求解过程，具体如下所示：

岩体在初始时刻处于平衡状态，且应力状态σ_ij及变形历史u_i已知；

施工开挖后岩体体积由Ω₀变为Ω，在应力边界S_σ上产生边界面力增量

在位移边界S_u上产生位移增量

使体积Ω内的应变增量Δε_ij、应力增量Δσ_ij和位移增量Δu_i满足下列控制方程及边界条件：

Δσ_ij,j+Δb_i＝0inΩ

其中，Δσ_ij,j为应力增量的偏导数，Δb_i为岩体体积力增量，

为弹性应变增量，

为非弹性应变增量，

为岩体的弹性刚度张量；n_j为应力边界S_σ的法向量；Δu_i,j、Δu_j,i为位移增量的偏导数；

采用有限单元法进行求解，设K₀、u₀、F₀分别为岩体开挖前初始应力场条件下的刚度矩阵、位移向量、开挖面节点力向量，下标0代表未开挖，则：

K₀u₀＝F₀

隧道开挖时，用K₁，u₁，F_exv1分别表示岩体开挖后的刚度矩阵、位移向量、开挖面节点力向量，下标1表示开挖，则：

K₁u₁＝F_exv1

其中，B^T为几何矩阵的对称矩阵；N^T为插值形函数矩阵的对称矩阵；b为边界面体力矢量；σ₀为岩体的初始应力；Ω_exv1为开挖掉的区域；

为边界面力；

岩体开挖后的开挖面节点力向量F_exv1即为在隧道开挖面节点上施加的等效开挖荷载。

步骤S4中，在开挖隧道江水渗流分析模型的基础上，采用TOUGHREACT软件进行初始渗流场分析，若洞周围岩中流体饱和且只有一种液相组分即水，则流体运动控制方程的表达形式为：

其中，ρ和μ分别为为水的密度和粘滞系数；g为重力加速度；q为水的单位体积质量源汇项；p为水压力；K为岩体渗透率；φ为孔隙率；t表示时间；

表示梯度算子；

采用TOUGHREACT软件进行渗流计算，得到计算区域中不同单元的中心点处水压力，根据不同单元中心点处的水压力，得到隧道开挖力学分析模型中各单元节点的水压力，具体如下所示：

隧道开挖力学分析模型中单元节点M的坐标为(x_M,y_M,z_M)，若与单元节点M相关联的单元有n个，根据TOUGHREACT软件的渗流计算结果，可知与单元节点M相关联的各单元中心点处的水压力为p_i，i＝1,2,…n，相关联的各单元中心点的坐标为(x_i,y_i,z_i)，i＝1,2,…n，则采用反距离插值法计算得到的隧道开挖力学分析模型中单元节点M的水压力p_M为：

步骤S4中，基于细观力学方法，得到水力耦合条件下饱和岩体的本构方程，具体如下所示：

对于饱和含任意微裂纹分布的岩体，在宏观应力Σ与水压力p共同作用下，基于Mori-Tanaka细观均匀化方法的自由焓表达式为：

其中，宏观应力Σ泛指岩体中应力大小，开挖荷载作用在隧道开挖面上，开挖荷载会引起宏观应力Σ的改变，进而引起岩体应力产生变化和发生变形；

p为隧道开挖力学分析模型中的水压力；

为微裂纹产生的宏观应变；

β表示微裂纹张开的内变量；

代表求两个向量并矢积对称分量的符合；

表示单位球面用来反映任意方向分布的微裂纹，n为微裂纹的单位法向量；

γ(n)表示单位法向量为n的微裂纹的滑移内变量矢量，γ表示微裂纹的滑移内变量；

S^s为各向同性固体基质的弹性柔度张量，由基质弹性模量E^s和泊松比v^s确定；

d表示岩体细观损伤的内变量，d＝Num×a³，Num为微裂纹密度表示单位体积的微裂纹数量，a为微裂纹的平均半径；

弹性常数H₀＝3E^s/{16[1-(v^s)²]}，H₁＝H₀(1-v^s/2)；

δ表示二阶单位张量；

B和N分别为Biot系数张量和Biot模量，B＝b₀δ，

b₀和φ₀分别为岩体各向同性Biot系数和初始孔隙率；

对岩体的自由焓关于宏观应力Σ与水压力p求导，可得水力耦合条件下饱和岩体的本构方程为：

其中，E表示宏观应变，φ为孔隙率；

洞周围岩内部微结构的变化分析如下所示：

水力耦合条件下饱和岩体的本构方程中，用于描述洞周围岩内部微结构的变化的内变量d、β和γ，分别根据对应的共轭热力学力F^d、F^β和F^γ确定，具体如下所示：

其中，热力学力F^β和F^γ分别代表微裂纹的局部法向和切向有效应力；

若微裂纹张开，即F^β＝0和F^γ＝0，则根据上式确定微裂纹张开的内变量β与滑移的内变量γ；

若微裂纹呈闭合状态，F^β＜0，则确定β和γ时考虑微裂纹的剪切滑移和法向压缩闭合效应，采用关联的Mohr-Coulomb准则，F＝|F^γ|+F^βtanφ_c，模拟微裂纹的滑移剪胀变形，法向闭合变形则利用双曲线模型表征，β＝-F^ββ₀/(k₀β₀-F^β)，其中，φ_c、β₀、k₀分别为微裂纹的内摩擦角、最大闭合量、初始法向刚度；

岩体细观损伤的内变量d，采用如下所示微裂纹损伤演化准则确定：

其中，V(d_c)为微裂纹损伤扩展抵抗力最大值，d_c为临界损伤，d_c与峰值应力对应的非弹性应变有关。

步骤S5中，根据洞周围岩内部微结构的变化，基于Voigt上限模型建立岩体渗透率演化方程，具体如下所示：

其中，K为岩体渗透率；φ为孔隙率；δ表示二阶单位张量；k^s为基质的各向同性渗透率；

为在初始损伤为d₀和体积率为β₀下的微裂纹的渗透率；参数χ用于表征微裂纹连通性随损伤演化的变化。

步骤S6中，采用交叉迭代算法进行水力耦合过程的求解，确定内变量d、β和γ，交叉迭代算法具体包括以下步骤：

S61，在每个时间步内，基于开挖隧道江水渗流分析模型，利用TOUGHREACT软件进行流体流动过程模拟，计算开挖隧道江水渗流分析模型中的水压力；

S62，将TOUGHREACT计算出的水压力代入隧道开挖力学分析模型中，进行考虑水压力变化影响的力学损伤过程的计算，分析洞周围岩内部微结构的变化，根据洞周围岩内部微结构的变化更新岩体渗透率K和孔隙率φ；

S63，将更新的岩体渗透率K和孔隙率φ代入TOUGHREACT软件中进行渗流计算，重新确定开挖隧道江水渗流分析模型中的水压力；

S64，重复上述步骤S62～S63，直到计算指定计算时长，得到江水渗流过程隧道开挖作用后围岩的应力Σ、应变E、位移u和水压力p，以及得到用于描述洞周围岩内部微结构的变化的内变量d、β和γ。

分别对不同江水位和不同岩体结构特性进行分析，得到不同江水位和不同岩体结构特性下江水渗流过程中的隧道应力变形。

本发明的优点在于：

(1)本发明方法利用基于细观力学的岩体水力耦合分析方法，对穿江隧道开挖过程中江水渗流的影响效应进行了***的分析，对隧道围岩细观结构特征进行了细致的模拟，对江水渗流过程与围岩损伤变形过程的耦合关系进行了合理设置，对穿江隧道开挖过程中江水渗流的影响效应进行了***的分析，应用前景良好。

(2)本发明方法开展江水渗流对穿江隧道开挖影响分析时，考虑隧道围岩微裂纹变化对围岩应变与渗透特性的改变，将岩体宏观力学响应与细观结构演化联系起来，模型参数的物理意义更加明确，具有创新意义，应用前景好。

(3)本发明还进一步分析了江水变化与岩体结构特征变化等因素在穿江隧道开挖施工中，对江底竖直变形的影响，具有创新意义，工程应用前景广泛。

附图说明

图1为本发明的江水渗流对穿江隧道开挖影响的分析方法流程图。

图2为本发明实施例的隧道初始有限元计算模型示意图。

图3为本发明实施例的隧道江水渗流有限元计算模型示意图。

图4为本发明实施例的反距离加权插值节点关联单元示意图。

图5为本发明实施例的岩体水力损伤耦合模拟示意图。

图6为本发明实施例的未考虑江水渗流时不同江水位下隧道开挖引起江底竖直变形。

图7为本发明实施例的考虑江水渗流时不同江水位下隧道开挖引起江底竖直变形。

图8为本发明实施例的不同微裂纹结构特征下隧道开挖引起江底竖直变形。

图9为本发明实施例的微裂纹各向同性时开挖隧道洞周水压力。

图10为本发明实施例的微裂纹优先垂直发育时开挖隧道洞周水压力。

图11为本发明实施例的微裂纹优先水平发育时开挖隧道洞周水压力。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

由图1所示，本发明提供江水渗流对穿江隧道开挖影响的分析方法，包括以下步骤：

S1，根据穿江隧道过江施工段，建立隧道开挖力学分析模型；其中，所述力学分析模型在进行初始应力设定时考虑不同江水位的影响。

建立隧道开挖力学分析模型，具体包括：

S11，建立初始有限元力学计算模型：

根据穿江隧道工程地质条件、隧道尺寸结构与布置以及江水条件，建立隧道初始有限元计算模型即初始模型，采用八节点六面体单元。模型顶部为江底，模型两侧与底部距隧道中心的尺寸大于5倍隧道洞径，整个模型均处于江水位以下，模型Z方向为竖直方向，围岩为饱和状态。

S12，根据初始模型边界方位，进行零位移法向约束边界条件的设定。

S13，根据实际隧道工程特点，采用初始自重应力计算或侧压力分析等方式，进行力学计算初始应力条件的设定。

初始自重应力计算方式：在模型顶部根据江水深度设置均布水压力荷载，然后进行模型在自重荷载作用下的应力计算，根据应力计算结果设置模型初始应力。

侧压力系数分析方式：根据侧压力系数K_c，采用如下计算公式设置竖直初始应力σ_z0与水平初始应力σ_x0＝σ_y0：

σ_z0＝γ_wH_w+γ_rH_r

σ_x0＝σ_y0＝K_cσ_z0

其中γ_w和γ_r为水和围岩的重度，H_w和H_r为江水深度和围岩的埋深。

S2，通过将隧道开挖区域挖除，建立隧道江水渗流有限元计算模型即江水渗流模型；

建立江水渗流模型，具体包括：

S21，在隧道初始有限元力学计算模型的基础上，通过将隧道开挖区域挖除，建立隧道有限元渗流计算模型。

S22，设置渗流分析边界与初始条件。

江水渗流模型的边界具体为：模型两侧和底部边界设为隔水边界，顶部取为固定水压力边界，其值根据江水位高程确定，开挖隧道洞周边界在隧道开挖后水压力降为1标准大气压。

江水渗流模型的初始条件具体为：初始时模型内部水压力从模型顶部往下按照静水压力呈梯形分布；其中，模型顶部水压力为：水密度×重力加速度×江水深度；模型底部水压力为：水密度×重力加速度×(江水深度+模型高度)；水压力从顶部到底部线性增加，分布规律同梯形。

S3，基于隧道开挖力学分析模型，进行隧道开挖过程的力学分析，得到隧道开挖面上的开挖荷载。

步骤S3中，为保持江水渗流与开挖作用下隧道变形与渗流双场交叉迭代过程有限元计算模型的一致性，通过在隧道开挖面节点上施加等效的开挖荷载来模拟隧道施工开挖过程；

隧道开挖面节点上所施加的开挖荷载的有限元求解过程，如下所示：

假设岩体在初始时刻处于平衡状态且应力状态σ_ij及变形历史u_i已知。施工开挖后岩体体积由Ω₀变为Ω，在应力边界S_σ上产生边界面力增量

在位移边界S_u上产生位移增量

三维情况边界都是面，作用力都施加在面上，称边界面力。

求解Ω内的应变增量Δε_ij、应力增量Δσ_ij和位移增量Δu_i，使之满足下列控制方程及边界条件：

Δσ_ij,j+Δb_i＝0inΩ

其中，Δσ_ij,j为应力增量的偏导数，Δb_i为岩体体积力增量，

为弹性应变增量，

为非弹性应变增量，

为岩体的弹性刚度张量；n_j为应力边界S_σ的法向量；Δu_i,j、Δu_j,i为位移增量的偏导数；

采用有限单元法对开挖问题进行求解时，设K₀、u₀、F₀分别为岩体开挖前初始应力场条件下的刚度矩阵、位移向量和等效节点力向量，下标0代表未开挖，则：

K₀u₀＝F₀

隧道开挖时，用K₁，u₁，F_exv1分别为岩体开挖后刚度矩阵、位移向量和等效节点力向量，下标1表示开挖第1步，则：

K₀u₀＝F₀

其中：B^T为几何矩阵的对称矩阵；N^T为插值形函数矩阵的对称矩阵；b为边界面体力矢量；σ₀为岩体的初始应力；Ω_exv1为开挖掉的区域；

为边界面力。

岩体开挖后的开挖面节点力向量F_exv1即为在隧道开挖面节点上施加的等效开挖荷载。

S4，在开挖隧道江水渗流分析模型的基础上，采用TOUGHREACT软件进行初始渗流场分析，获得计算区域不同单元中心点的水压力大小，然后通过插值分析由单元中心点水压力结果得到隧道开挖力学分析模型中各单元节点的水压力值；采用基于细观力学的损伤力学模型，开展水压力和开挖荷载作用下的隧道应力变形分析，得到洞周围岩内部微结构的变化；

洞周围岩指：开挖隧道周围的岩体。内部微结构指：岩体内部微小结构面，对岩体变形破坏与渗流过程起控制作用。

步骤S4中，利用TOUGHREACT软件进行渗流计算，其流体运动控制方程一般表达形式为：

其中，

代入上式后，得到：

其中，M_κ为单位体积组分κ的质量；S_l和ρ_l分别为流体饱和度和密度，

为液相组分κ的质量分数；q_κ为组分κ的单位体积质量源汇项；F_κ为组分κ的质量通量；K为岩体渗透率；φ为孔隙率；p_l表示液相流体压力；div表示散度算子；

表示梯度算子；k_rl表示液相的相对渗透系数；下标r、l分别表示相对渗透系数与液相；t表示时间。

当洞周围岩中流体饱和且只有一种液相组分，如水，k_rl＝1、

则上式流体运动控制方程可写成：

其中，K为岩体渗透率；ρ和μ为水密度和粘滞系数；g为重力加速度；q为单位体积质量源汇项；p为水压力；φ为孔隙率。

TOUGHREACT渗流计算得到的水压力计算结果位于单元中心点，而有限元力学过程需用到单元节点处的水压力，在进行力学过程求解时，通过反距离加权平均法将单元中心的水压力值插值到单元节点上。若与节点M相关联的单元有n个，各单元中心的水压力值为p_i,i＝1,2,…n，节点M的坐标为(x_M,y_M,z_M)，关联单元中心的坐标为(x_i,y_i,z_i)，i＝1,2,…n，则采用反距离插值法计算得到的节点M的水压力p_M为

数值模拟采用的水力耦合条件下岩体本构模型考虑岩体细观结构特性，基于细观力学方法获得。

对饱和含任意微裂纹分布的岩体，其在宏观应力Σ与水压力p共同作用下基于Mori-Tanaka细观均匀化方法的自由焓表达式为：

其中，宏观应力Σ泛指岩体中应力大小。开挖荷载只作用在局部区域，开挖荷载会引起宏观应力的改变，进而岩体发生应变、变形。

为微裂纹产生的宏观应变；

β代表微裂纹张开的内变量；

代表求两个向量并矢积对称分量的符合；

表示单位球面用来反映任意方向分布的微裂纹，n为微裂纹的单位法向量；

γ(n)表示单位法向量为n的微裂纹的滑移内变量矢量；γ表示微裂纹的滑移内变量；

S^s为各向同性固体基质的弹性柔度张量，由基质弹性模量E^s和泊松比v^s确定；

d＝Num×a³表示岩体细观损伤的内变量，Num为微裂纹密度表示单位体积的微裂纹数量，a为微裂纹的平均半径；

弹性常数H₀＝3E^s/{16[1-(v^s)²]}，H₁＝H₀(1-v^s/2)；

δ表示二阶单位张量；

B和N分别为Biot系数张量和Biot模量，B＝b₀δ，

b₀和φ₀分别为岩体各向同性Biot系数和初始孔隙率。

宏观应力Σ泛指岩体中应力大小。开挖荷载只作用在局部区域，开挖荷载会引起宏观应力的改变，进而岩体发生应变、变形。

对岩体自由焓关于宏观应力Σ和p求导，可得水力耦合条件下饱和岩体的本构方程为：

其中，E表示宏观应变，φ为孔隙率；

本构方程中描述岩体内部微裂纹细观结构变化的内变量d、β和γ的演化根据其对应的共轭热力学力F^d、F^β和F^γ确定

其中，热力学力F^β和F^γ实质上代表微裂纹的局部法向和切向有效应力。当微裂纹张开时，有F^β＝0和F^γ＝0，此时可根据上式确定β和γ，当微裂纹呈闭合状态时(F^β＜0)，确定β和γ时考虑微裂纹的剪切滑移和法向压缩闭合效应，采用关联的Mohr-Coulomb准则，F＝|F^γ|+F^βtanφ_c，模拟微裂纹的滑移剪胀变形，而法向闭合变形则利用双曲线模型表征，β＝-F^ββ₀/(k₀β₀-F^β)，其中φ_c、β₀和k₀分别为微裂纹的内摩擦角、最大闭合量和初始法向刚度。内变量d的演化采用如下的微裂纹损伤演化准则：

式中，V(d_c)为微裂纹损伤扩展抵抗力最大值，d_c为临界损伤与峰值应力对应的非弹性应变有关。

S5，根据岩体内部微结构变化，采用基于细观均匀化方法得到的岩体渗透性演化方程进行洞周围岩渗透张量的计算，并将其传递到渗流分析模型中。

步骤S5中，数值模拟考虑岩体渗透特性变化，且其与微裂纹细观结构变化密切相关，采用基于Voigt上限模型建立的渗透张量演化方程进行估计：

其中，k^s为基质的各向同性渗透率；

为初始损伤和体积率分别为d₀和β₀时微裂纹的渗透率；参数χ用于表征微裂纹连通性随损伤演化的变化，其值可通过岩体渗透性演化试验数据反算确定。

S6，采用双场交叉迭代算法进行江水渗流过程开挖隧道围岩应力、变形与渗流过程的求解。

步骤S6中，采用双场交叉迭代算法进行耦合过程的求解，主要包括流体流动计算模块和力学损伤耦合计算模块，其中流体流动计算采用广泛应用的TOUGHREACT多场多相流体模拟程序，力学损伤耦合计算基于预测-校正算法采用增量形式的有限元代码进行力学本构方程求解，然后确定各表征微裂纹细观结构的内变量d、β和γ。迭代算法具体包括：

S61，在每个时间步内，先采用TOUGHREACT进行流体流动过程模拟，确定开挖隧道江水渗流分析模型中的水压力分布结果。

S62，将TOUGHREACT渗流计算出的开挖隧道江水渗流分析模型中的水压力代入隧道开挖力学分析模型中，进行考虑水压力变化影响的力学损伤过程的计算，分析洞周围岩内部微结构的变化，根据洞周围岩内部微结构的变化更新岩体渗透率K和孔隙率φ。

S63，将更新的岩体渗透率K和孔隙率φ代入TOUGHREACT软件中进行渗流计算，重新确定开挖隧道江水渗流分析模型中的水压力分布结果。

S64，重复上述步骤S62～S63，直到计算指定计算时长，得到江水渗流过程隧道开挖作用后围岩的应力Σ、应变E、位移u和水压力p，以及用于描述洞周围岩内部微结构的变化的内变量d、β和γ，据此对穿江隧道的开挖扰动效应进行评价分析，为工程设计施工提供支撑。

S7、考虑江水位和岩体结构特性变化情况开展数值分析，得到不同江水位和岩体结构特性下江水渗流过程中的隧道应力变形，进一步得到江水位和岩体结构对江底沉降的影响情况；

步骤S7中，可根据江水位变化范围与岩体结构特征，模拟不同江水位和微裂纹结构特性情况下隧道开挖引起的江底竖直变形。

实施例：

本实例选取安徽省首条过江隧道——芜湖超大直径过江隧道穿越长江江心段，借助自主开发的基于细观力学的岩体水力耦合数值模型，开展江底隧道开挖过程的水力耦合响应模拟。

如图1所示，本实例的江水渗流对穿江隧道开挖影响的分析方法，包括以下步骤：

S1，建立隧道开挖力学分析模型，具体如下所示：

S11，根据隧道几何模型，建立隧道有限元计算模型：

根据芜湖过江隧道穿越长江江心段的布置特点，建立隧道初始有限元计算模型即初始模型，如图2所示，模型六面体网格数为7207，节点数为14600，模型平面尺度为200m×200m，模型顶部为江底，其高程为-32.8m，整个模型均处于长江水位以下，隧道顶部距江底13.57m，围岩为饱和状态。

S12，施加位移边界条件：

对隧道初始有限元计算模型的底部和两侧均施加零位移法向约束，如图2所示。

S13，给定初始应力条件：

根据长江水头条件，采用侧压力系数分析方式给定初始应力条件：竖直初始应力σ_z0＝γ_wH_w+γ_rH_r，水平初始应力σ_x0＝σ_y0＝K_cσ_z0，其中，γ_w和γ_r为水和围岩的重度，H_w和H_r为江水深度和围岩的埋深，K_c为侧压力系数，K_c取值为0.5。

S2，建立建立隧道江水渗流有限元计算模型即江水渗流模型，具体如下所示：

S21，建立隧道江水渗流有限元计算模型：

根据芜湖过江隧道穿越长江江心段的布置特点，建立隧道江水渗流有限元计算模型，如图3所示，模型六面体网格数为5406，节点数为11120。模型平面尺度为200m×200m，模型顶部为江底，其高程为-32.8m，整个模型均处于长江水位以下，隧道顶部距江底13.57m，围岩为饱和状态。

S22，设置渗流分析边界与初始条件：

由图3所示，计算模型两侧设为隔水边界，顶部取为固定水压力边界，其值根据长江水位高程确定，如江水位为4m时，等效压力为0.368MPa，初始时模型内部水压力呈梯形分布，且水力耦合模拟过程中隧道开挖后，洞周围岩水压力降为1标准大气压(0.1MPa)，水力耦合计算时长为1d。

S3，基于隧道开挖力学分析模型，进行隧道开挖过程的力学分析，得到隧道开挖面上的开挖荷载；为保持江水渗流与开挖作用下隧道变形与渗流双场交叉迭代过程有限元计算模型的一致性，通过在隧道开挖面节点上施加等效的开挖荷载来模拟隧道施工开挖过程。

S4，在开挖隧道江水渗流分析模型的基础上，采用TOUGHREACT软件进行初始渗流场分析，获得计算区域不同单元中心点的水压力大小，然后通过插值分析由单元中心点水压力结果得到隧道开挖力学分析模型中各单元节点的水压力值；采用基于细观力学的损伤力学模型，开展水压力和开挖荷载作用下的隧道应力变形分析，得到洞周围岩内部微结构的变化。

在开挖隧道江水渗流分析模型的基础上，采用TOUGHREACT软件进行初始渗流场分析，若洞周围岩中流体饱和且只有一种液相组分即水，则流体运动控制方程的表达形式为：

其中，ρ和μ为水的密度和粘滞系数；g为重力加速度；q为水的单位体积质量源汇项；p为水压力；K为岩体渗透率；φ为孔隙率；t表示时间；

表示梯度算子；

采用TOUGHREACT软件进行渗流计算，得到计算区域中不同单元的中心点处水压力，根据不同单元中心点处的水压力，得到隧道开挖力学分析模型中各单元节点的水压力，具体如下所示：

如图4所示，隧道开挖力学分析模型中单元节点M的坐标为(x_M,y_M,z_M)，若与单元节点M相关联的单元有n个，根据TOUGHREACT软件的渗流计算结果，可知与单元节点M相关联的各单元中心点处的水压力为p_i，i＝1,2,…n，相关联的各单元中心点的坐标为(x_i,y_i,z_i)，i＝1,2,…n，则采用反距离插值法计算得到的隧道开挖力学分析模型中单元节点M的水压力p_M为：

S5，根据岩体内部微结构变化，采用基于细观均匀化方法得到的岩体渗透性演化方程进行洞周围岩渗透张量的计算，并将其传递到开挖隧道江水渗流分析模型中。

S6，将岩体渗透率代入开挖隧道江水渗流分析模型中，继续进行渗流场分析，更新开挖隧道江水渗流分析模型中的水压力，如图5所示，按照步骤S4～S5的方式，采用交叉迭代算法，进行江水渗流过程开挖隧道围岩应力、变形与渗流过程的求解，得到江水渗流过程隧道开挖作用后围岩的应力Σ、应变E、位移u和水压力p，以及用于描述洞周围岩内部微结构的变化的内变量d、β和γ，据此对穿江隧道的开挖扰动效应进行评价分析，为工程设计施工提供支撑。

S7，考虑江水位和岩体结构特性变化情况开展数值分析，得到不同江水位和岩体结构特性下的江底地表沉降数值，进一步得到江水位和围岩结构对江底沉降的影响情况。

步骤S7的数值分析具体如下所示：

S71，岩体物理力学参数：

根据现场勘测和室内试验成果，并参考相关经验进行选取穿江隧道围岩计算参数如下：基质弹性模量E^s＝500MPa、基质泊松比v^s＝0.25、损伤扩展抵抗力V(d_c)＝1×10^-3MPa、损伤变量临界值d_c＝1.0、岩块Biot系数b₀＝0.8、岩体重度γ＝25KN/m³、基质渗透率k^s＝×10^-16m²、微裂纹初始渗透率k₀ ^c＝1.5×10^-12m²。

S72，数值模拟结果

以本实例的隧道工程为例，模拟当长江水位分别为0m、4m、8m和12m四种情况下，隧道开挖对地表竖直变形的影响。此时，岩体内部初始微裂纹假定呈各向同性分布。

在未考虑与考虑江水渗流作用下不同江水位下隧道开挖引起江底竖直变形，如图6-7所示。由图可知，不考虑江水渗流作用时，隧道开挖荷载引起的最大竖直变形位于隧道中轴线区域，江底地面产生沉降，远离隧道区域地面发生少量***，江水位分别为0m、4m、8m和12m时，江底地表沉降最大值分别为0.0067m、0.0076m、0.0085m、0.0094m，地表***最大值分别为0.0045m、0.0044m、0.0043m、0.0042m，江水位变化引起的沉降增幅不超过0.003m。考虑江水渗流作用后，江底地面沉降明显，无***现象，隧道中轴线处江底地面沉降最大，在0m江水位时为0.036m，12m江水位时为0.049m，江水变化引起的增幅达0.013m。由此可见，江水渗流作用对江底地面沉降变形影响很大，且随着江水位增大，影响效应越明显，实际分析中考虑江水渗流作用非常重要。

岩体力学特性和渗透性受内部初始微裂纹分布影响，受长期地质作用影响，岩体内部微裂纹结构复杂。本实例模拟分析岩体内部微裂纹呈各向异性分布对隧道开挖引起的江底地面竖直变形，考虑微裂纹优先呈水平与垂直发育两种工况，江水位取4m。

不同初始微裂纹分布下隧道开挖引起江底竖直变形如图8所示。由图可知，初始微裂纹分布特征对江底隧道开挖引起的地表变形响应影响较大，微裂纹呈各向同性分布时隧道中轴线处地表沉降值最大，优先水平发育时次之，垂直发育时最小。微裂纹优先垂直发育时，隧道开挖引起江底地表竖直沉降变形呈抛物线分布，与各向同性情况时规律一致，但沉降值变小，最大沉降值为0.0327m。与初始微裂纹呈各向同性分布的结果对比，优先水平发育时，江底地表竖直沉降最大值减小，为0.0372m，但远离隧道中轴线2倍洞直径区域，竖直沉降增大，隧道开挖引起的江底地表沉降范围扩大。

渗流对开挖隧道围岩变形的影响与水压力分布密切相关。不同初始微裂纹分布下开挖隧道洞周水压力如图9-11所示。由图可知，微裂纹优先水平发育时，水压力降低区域明显增大，从而导致远离开挖隧道较远区域地表也产生较大的沉降变形。由此可见，初始微裂纹各向异性分布对江底隧道开挖引起的地表变形与水压力分布影响显著，岩体内部微裂纹等细观结构特征的准确表征对江底隧道开挖水力耦合响应模拟非常重要。

因此，本发明通过分析江水渗流对穿江隧道开挖的影响，以便合理进行江底隧道开挖施工保护，应用前景良好。

以上仅为本发明创造的较佳实施例而已，并不用以限制本发明创造，凡在本发明创造的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明创造的保护范围之内。

除非特别说明，本发明文本中的以前后顺序出现的各个流程并不必然存在先后的执行顺序。尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (9)

1.江水渗流对穿江隧道开挖影响的分析方法，其特征在于，包括以下步骤：

S1，建立隧道开挖力学分析模型；

S2，基于隧道开挖力学分析模型，通过将隧道开挖区域挖除，建立开挖隧道江水渗流分析模型；

S3，基于隧道开挖力学分析模型，进行隧道开挖过程的力学分析，得到隧道开挖面上的开挖荷载；

S4，基于开挖隧道江水渗流分析模型，进行初始渗流场分析，分析开挖隧道江水渗流分析模型中的水压力；将开挖隧道江水渗流分析模型中水压力的分析结果代入隧道开挖力学分析模型中，基于隧道开挖力学分析模型，进行开挖载荷和水压力共同作用下的隧道应力与变形分析，分析洞周围岩内部微结构的变化；

S5，建立岩体渗透率演化方程，并根据洞周围岩内部微结构的变化，计算洞周围岩的岩体渗透率；

S6，将岩体渗透率代入开挖隧道江水渗流分析模型中，继续进行渗流场分析，更新开挖隧道江水渗流分析模型中的水压力，按照步骤S4～S5的方式，采用交叉迭代算法，分析江水渗流过程中隧道的应力与变形，以及洞周围岩内部微结构的变化。

2.根据权利要求1所述的江水渗流对穿江隧道开挖影响的分析方法，其特征在于，步骤S1中，建立隧道开挖力学分析模型，具体包括以下步骤：

S11，根据穿江隧道工程地质条件、隧道尺寸结构与布置以及江水条件，采用八节点六面体单元，建立隧道初始有限元计算模型即初始模型；其中，初始模型的顶部为江底，整个模型均处于江水位以下，Z方向为竖直方向，围岩为饱和状态；

S12，对初始模型的底部和两侧均施加零位移法向约束；

S13，根据实际隧道工程特点，采用初始自重应力计算或采用侧压力系数分析的方式，设置初始模型的初始应力；

所述初始自重应力计算的具体方式为：根据江水深度在初始模型顶部设置均匀分布的水压力荷载，在隧道自重荷载作用下对初始模型的应力进行计算，根据应力计算结果设置初始模型的初始应力；

所述侧压力系数分析的具体方式为：根据侧压力系数K_c，采用如下计算公式，对初始模型的的初始应力进行设置：

σ_z0＝γ_wH_w+γ_rH_r

σ_x0＝σ_y0＝K_cσ_z0

其中，σ_z0为竖直方向的初始应力，σ_x0为水平面沿x方向的初始应力，σ_y0为水平面沿y方向的初始应力，γ_w和γ_r分别为水和围岩的重度，H_w和H_r分别为江水深度和围岩埋深。

3.根据权利要求1所述的江水渗流对穿江隧道开挖影响的分析方法，其特征在于，步骤S2中，建立开挖隧道江水渗流分析模型，具体包括以下步骤：

S21，在隧道开挖力学分析模型的基础上，将隧道开挖区域挖除，建立隧道江水渗流有限元计算模型即江水渗流模型；

S22，设置江水渗流模型的边界与初始条件；

江水渗流模型的边界具体为：模型两侧和底部边界设为隔水边界；顶部设为固定水压力边界，顶部的固定水压力边界值根据江水位高程确定；在隧道开挖后，隧道洞周边界水压力降为1标准大气压；

江水渗流模型的初始条件具体为：初始时，江水渗流模型的内部水压力从顶部到底部呈梯形分布。

4.根据权利要求1所述的江水渗流对穿江隧道开挖影响的分析方法，其特征在于，步骤S3中，通过在隧道开挖面节点上施加等效的开挖荷载来模拟隧道施工开挖过程；

隧道开挖面节点上所施加的开挖荷载的有限元求解过程，具体如下所示：

岩体在初始时刻处于平衡状态，且应力状态σ_ij及变形历史u_i已知；

施工开挖后岩体体积由Ω₀变为Ω，在应力边界S_σ上产生边界面力增量

在位移边界S_u上产生位移增量

使体积Ω内的应变增量Δε_ij、应力增量Δσ_ij和位移增量Δu_i满足下列控制方程及边界条件：

Δσ_ij,j+Δb_i＝0inΩ

其中，Δσ_ij,j为应力增量的偏导数，Δb_i为岩体体积力增量，

为弹性应变增量，

为非弹性应变增量，

为岩体的弹性刚度张量；n_j为应力边界S_σ的法向量；Δu_i,j、Δu_j,i为位移增量的偏导数；

采用有限单元法进行求解，设K₀、u₀、F₀分别为岩体开挖前初始应力场条件下的刚度矩阵、位移向量、开挖面节点力向量，下标0代表未开挖，则：

K₀u₀＝F₀

隧道开挖时，用K₁，u₁，F_exv1分别表示岩体开挖后的刚度矩阵、位移向量、开挖面节点力向量，下标1表示开挖，则：

K₁u₁＝F_exv1

其中，B^T为几何矩阵的对称矩阵；N^T为插值形函数矩阵的对称矩阵；b为边界面体力矢量；σ₀为岩体的初始应力；Ω_exv1为开挖掉的区域；

为边界面力；

岩体开挖后的开挖面节点力向量F_exv1即为在隧道开挖面节点上施加的等效开挖荷载。

5.根据权利要求1所述的江水渗流对穿江隧道开挖影响的分析方法，其特征在于，步骤S4中，在开挖隧道江水渗流分析模型的基础上，采用TOUGHREACT软件进行初始渗流场分析，若洞周围岩中流体饱和且只有一种液相组分即水，则流体运动控制方程的表达形式为：

其中，ρ和μ分别为为水的密度和粘滞系数；g为重力加速度；q为水的单位体积质量源汇项；p为水压力；K为岩体渗透率；φ为孔隙率；t表示时间；

表示梯度算子；

采用TOUGHREACT软件进行渗流计算，得到计算区域中不同单元的中心点处水压力，根据不同单元中心点处的水压力，得到隧道开挖力学分析模型中各单元节点的水压力，具体如下所示：

隧道开挖力学分析模型中单元节点M的坐标为(x_M,y_M,z_M)，若与单元节点M相关联的单元有n个，根据TOUGHREACT软件的渗流计算结果，可知与单元节点M相关联的各单元中心点处的水压力为p_i，i＝1,2,…n，相关联的各单元中心点的坐标为(x_i,y_i,z_i)，i＝1,2,…n，则采用反距离插值法计算得到的隧道开挖力学分析模型中单元节点M的水压力p_M为：

6.根据权利要求1所述的江水渗流对穿江隧道开挖影响的分析方法，其特征在于，步骤S4中，基于细观力学方法，得到水力耦合条件下饱和岩体的本构方程，具体如下所示：

对于饱和含任意微裂纹分布的岩体，在宏观应力Σ与水压力p共同作用下，基于Mori-Tanaka细观均匀化方法的自由焓表达式为：

其中，宏观应力Σ泛指岩体中应力大小，开挖荷载作用在隧道开挖面上，开挖荷载会引起宏观应力Σ的改变，进而引起岩体应力产生变化和发生变形；

p为隧道开挖力学分析模型中的水压力；

为微裂纹产生的宏观应变；

β表示微裂纹张开的内变量；

代表求两个向量并矢积对称分量的符合；

表示单位球面用来反映任意方向分布的微裂纹，n为微裂纹的单位法向量；

γ(n)表示单位法向量为n的微裂纹的滑移内变量矢量，γ表示微裂纹的滑移内变量；

S^s为各向同性固体基质的弹性柔度张量，由基质弹性模量E^s和泊松比v^s确定；

d表示岩体细观损伤的内变量，d＝Num×a³，Num为微裂纹密度表示单位体积的微裂纹数量，a为微裂纹的平均半径；

弹性常数H₀＝3E^s/{16[1-(v^s)²]}，H₁＝H₀(1-v^s/2)；

δ表示二阶单位张量；

B和N分别为Biot系数张量和Biot模量，B＝b₀δ，

b₀和φ₀分别为岩体各向同性Biot系数和初始孔隙率；

对岩体的自由焓关于宏观应力Σ与水压力p求导，可得水力耦合条件下饱和岩体的本构方程为：

其中，E表示宏观应变，φ为孔隙率；

洞周围岩内部微结构的变化分析如下所示：

水力耦合条件下饱和岩体的本构方程中，用于描述洞周围岩内部微结构的变化的内变量d、β和γ，分别根据对应的共轭热力学力F^d、F^β和F^γ确定，具体如下所示：

其中，热力学力F^β和F^γ分别代表微裂纹的局部法向和切向有效应力；

若微裂纹张开，即F^β＝0和F^γ＝0，则根据上式确定微裂纹张开的内变量β与滑移的内变量γ；

若微裂纹呈闭合状态，F^β＜0，则确定β和γ时考虑微裂纹的剪切滑移和法向压缩闭合效应，采用关联的Mohr-Coulomb准则，F＝|F^γ|+F^βtanφ_c，模拟微裂纹的滑移剪胀变形，法向闭合变形则利用双曲线模型表征，β＝-F^ββ₀/(k₀β₀-F^β)，其中，φ_c、β₀、k₀分别为微裂纹的内摩擦角、最大闭合量、初始法向刚度；

岩体细观损伤的内变量d，采用如下所示微裂纹损伤演化准则确定：

其中，V(d_c)为微裂纹损伤扩展抵抗力最大值，d_c为临界损伤，d_c与峰值应力对应的非弹性应变有关。

7.根据权利要求6所述的江水渗流对穿江隧道开挖影响的分析方法，其特征在于，步骤S5中，根据洞周围岩内部微结构的变化，基于Voigt上限模型建立岩体渗透率演化方程，具体如下所示：

其中，K为岩体渗透率；φ为孔隙率；δ表示二阶单位张量；k^s为基质的各向同性渗透率；

为在初始损伤为d₀和体积率为β₀下的微裂纹的渗透率；参数χ用于表征微裂纹连通性随损伤演化的变化。

8.根据权利要求7所述的江水渗流对穿江隧道开挖影响的分析方法，其特征在于，步骤S6中，采用交叉迭代算法进行水力耦合过程的求解，确定内变量d、β和γ，交叉迭代算法具体包括以下步骤：

S61，在每个时间步内，基于开挖隧道江水渗流分析模型，利用TOUGHREACT软件进行流体流动过程模拟，计算开挖隧道江水渗流分析模型中的水压力；

S62，将TOUGHREACT计算出的水压力代入隧道开挖力学分析模型中，进行考虑水压力变化影响的力学损伤过程的计算，分析洞周围岩内部微结构的变化，根据洞周围岩内部微结构的变化更新岩体渗透率K和孔隙率φ；

S63，将更新的岩体渗透率K和孔隙率φ代入TOUGHREACT软件中进行渗流计算，重新确定开挖隧道江水渗流分析模型中的水压力；

S64，重复上述步骤S62～S63，直到计算指定计算时长，得到江水渗流过程隧道开挖作用后围岩的应力Σ、应变E、位移u和水压力p，以及得到用于描述洞周围岩内部微结构的变化的内变量d、β和γ。

9.根据权利要求1-8任意一项所述的江水渗流对穿江隧道开挖影响的分析方法，其特征在于，分别对不同江水位和不同岩体结构特性进行分析，得到不同江水位和不同岩体结构特性下江水渗流过程中的隧道应力变形。

Images