水工隧洞围岩与衬砌相互作用计算方法的比较及应用研究：基于摩擦滑移接触的统一分析框架与工程优化策略

本文系统建立了水工隧洞围岩-衬砌相互作用的统一分析框架，集成粘结、光滑和摩擦滑移三种接触模型，通过解析法与数值法对比揭示了接触行为机理。研究发现摩擦系数阈值（fr=0.59）控制接触模式转变，界面滑移使衬砌拉应力降低9.5%；弹性模量比E2/E1是应力敏感性的主导因素（2增至5时拱顶拉力提升80.3%）；侧压力系数λ与位移释放系数η显著影响应力分布模式。提出的变量降维算法将计算速度提升75%、精度提高2.2%，为水电工程安全性与经济性提供了量化优化方案。

　　

1 引言

中国西南山区峡谷地区通过流域梯级开发广泛利用丰富的水电资源。随着水电工程向高海拔地区转移，引水隧洞呈现规模更大、埋深更大的特点，对结构提出更复杂的需求。与交通隧道不同，水工隧洞承受内外水压力联合作用，使其在多场荷载下的结构完整性成为关键问题。根据基于连续介质力学的新奥法（NATM），隧道稳定性根本上通过衬砌与围岩的协同相互作用来维持。

围岩-衬砌界面相互作用模型的选择显著影响结构计算。尽管解析和数值方法在隧道工程中得到广泛应用，但三个持续存在的局限性阻碍了实际应用：衬砌-岩石界面材料非均质性导致的接触复杂性、推导快速解析解的计算效率不足，以及参数敏感性（如相对衬砌刚度E₂/E₁、厚度d/R₀、侧压力系数λ和位移释放系数η）量化不足。这些差距在摩擦滑动接触条件下尤为明显，成熟的 methodologies 仍然不发达。

先前研究主要沿着三种接触模型范式进展。圆形衬砌隧洞的分析通常采用弹性力学中的厚壁圆筒理论。这种方法主要区分两类问题：完全粘结接触和光滑接触问题。对于完全粘结接触，以完全位移连续性为特征，Bobet首次推导了考虑孔隙水效应的压力隧洞闭式解，后来扩展到浅埋隧洞。虽然Li等人研究了初始应力场下的弹性响应，但忽略了支护延迟——这一遗漏后来通过位移释放系数得到解决。考虑孔隙水压力的影响，Guo等人使用Mohr-Coulomb准则评估了全接触条件下衬砌隧洞围岩的稳定性，但他们仍然假设衬砌在开挖后立即建造。Wang等人考虑了支护延迟，但仅在径向，没有解决切向位移的影响。Han等人进一步将支护滞后集成到复变公式中。

相比之下，光滑接触模型施加零界面剪应力，最初使用相对刚度方法建立，并由Moore和Booker建立了严格的基准。最近的创新包括Gao的双层圆筒解和Lu的结合支护延迟与纯滑移条件的模型。

关于摩擦滑动接触，刚体的库仑摩擦模型简单、易用且应用广泛。Meschke等人采用库仑摩擦模型模拟桩与土的相互作用。分析复杂性源于非线性边界约束，导致大多数研究转向数值近似，如增广拉格朗日方法或分段弹簧系统。尽管Lu等人推导了圆形隧洞的解析解，但未进行所有三种接触模式的比较分析——特别是由不同工程参数调制的应力转换量化。Ahn和Pouya建立了衬砌与地层之间弹性接触的新解析解，并与其他工作进行了各种滑移条件的比较，但衬砌被建模为弹性圆形壳（二维中的梁）。

对衬砌隧洞现有研究的回顾揭示，最常见的方法假设围岩与衬砌之间完全粘结或完全光滑界面。然而，这些研究常常忽略开挖过程或未能考虑支护延迟的实际现象。在摩擦滑移接触条件下，Ahn和Pouya的工作将衬砌简化为梁，并未考虑衬砌安装前的开挖荷载比，同样假设衬砌的应用是瞬时的。虽然Lu考虑了支护延迟，但其优化模型的建立和求解过程存在固有缺陷，导致该方法难以扩展到任意截面形状的隧洞。因此，两个关键的研究缺陷仍然存在：缺乏比较粘结、光滑和摩擦滑移接触应力分布的统一框架，以及由于未量化的参数影响而导致的有限工程适用性。

2 岩-衬相互作用的解析方法

2.1 基本理论

圆形隧洞中围岩和衬砌的应力应变计算由于其高长径比和高对称性特点，通常简化为关键截面求解。围岩与衬砌之间的接触可分为三种类型：完全粘结接触、光滑接触和摩擦滑移接触。

完全粘结接触假设界面上任意点无相对位移。在数值方法中，这通常通过共享节点模拟，假设法向应力、切向剪应力和位移的连续性（下标I和II分别代表围岩和衬砌）。应力和位移连续性条件分别由方程(1)和(2)表示。

光滑接触保持法向应力和位移连续性，但允许切向位移不连续且剪应力为零。此时，边界条件和连续性条件的表达式由方程(3)和(4)给出。

摩擦滑移接触代表一种中间状态，由库仑摩擦理论控制，如方程(5)所示。当摩擦系数为零或超过某个阈值时，可以推导出两种极端接触条件（完全结合接触和光滑接触）的解析解。这种方法也常被称为基于复变函数理论的优化算法，其中，当摩擦系数低于临界阈值时可能发生相对滑动，而较高值则诱发粘结行为。因此，应用复变函数理论需要根据边界条件和接触模型确定解析函数系数，以推导应力和变形场。

2.2 解析函数的工程构建

建立了深埋圆形隧洞在初始应力和内水压力下的计算模型。参数指定为：p₀ = 2.0 MPa, R₁ = 3.0 m, R₀ = 2.5 m, μ₁ = 0.25 (岩石), μ₂ = 0.20 (衬砌), E₂/E₁ = 10.0, η = 0.20。系数η量化支护延迟，其中η = 0表示开挖后立即安装衬砌，η = 1对应于岩体完全变形后安装。对于摩擦滑移接触，初始采用f_r = 0.5。

隧洞开挖后，边界上的面力被释放，由此开挖产生的位移可以用方程(6)表示。

围岩和衬砌中的两个复势函数是γ_R(z), χ_R(z), γ_L(z), χ_L(z)，它们可以用Taylor和Laurent级数表示。围岩的应力和位移表达式由方程(7)和(8)给出。

在衬砌中，衬砌的应力和位移表达式由方程(9)和(10)给出。

基于2.1节描述的边界条件和接触条件，可以制定一个方程组（或不等式）来确定解析函数的系数。三种不同接触模式下解析解的计算过程分别如下：

1. 1.
   完全粘结接触需要考虑衬砌的内部边界条件、接触界面处的连续性条件、位移边界条件以及接触界面条件。对于圆形隧洞，根据推导过程，仅需要九个参数来确定相应的解析函数。
2. 2.
   光滑接触修改了五个方程，同时由于接触界面连续性和接触条件的差异，保留了方程(11)-(13)和(17)。
3. 3.
   摩擦滑移接触由于缺乏接触条件方程，只能列出八个方程。并且与光滑接触不同的方程由方程(25)给出。

有九个未知数但只有八个方程，因此需要使用附加条件进行最优求解。

摩擦滑动接触条件的特征是通过识别最小相对滑动的发生点作为临界点。根据库仑摩擦理论，相应的不等式方程如下。

以接触面上相对滑动最小化为目标函数。

最后，通过采用优化计算方法，也可以获得满足边界条件的解析函数系数。

这里，考虑了一种新的优化方法来解决这个问题。不等式约束和目标函数的表达式分别在方程(26)和(27)中给出。原始优化模型包含9个设计变量和8组等式方程。通过选择x_i (i = 1, … , 9)中的任何一个变量作为设计变量——例如，x₉——并为其分配一个初始值，剩余的8个变量x_i (i = 1, … , 8)可以通过求解8个方程来确定。随后，将这些x_i (i = 1, … , 8)的值代入以下优化模型以获得优化结果。

类似地，通过采用增广罚函数法进行优化，并且与包含等式约束的传统罚函数形式相比，方程(28)中给出的罚函数形式不包含带有∑_j=1⁸a_j²(X)的二次损失项。而且，初始点仅涉及x₉，这大大简化了优化的数学模型。

为了比较两种模型的准确性，我们以摩擦滑移接触为例。当f_r = 0.1时，新旧模型计算的目标函数值（即切向位移差）分别为2.4129和2.4659，提高了约2.2%。比较两种模型获得的应力结果，并以衬砌内边界处的径向应力为例，理论值为2.0 MPa。使用新优化模型获得的应力的最大相对误差为2.025 × 10^?2%，而使用旧优化模型获得的最大相对误差为2.047%。因此，相对于理论解，旧优化模型的结果表现出较大的误差。同时，在旧优化模型中，目标函数被访问约1000至3000次，而在新模型中，这个数字已减少到200-300的范围。这种改进归因于本文提出的方法可以通过方程求解确保满足a_j(X)=0，使其比以前的优化模型更准确。此外，设计变量的数量大大减少，显著提高了优化速度。而且，该方法可以扩展到解决任意截面隧洞的优化问题。对于复杂截面，未知数的数量可达6N + 4，其中N通常大于60。当使用原始优化模型时，求解过程常常陷入局部最优甚至无法收敛。通过应用新模型，未知数的系数可以减少约80%，有效防止了过多未知数导致问题无法求解的问题。

2.3 求解解析函数的过程

确定了解析函数的系数后，衬砌内的应力分量（σ_r^L; σ_θ^L, τ_rθ^L)由方程(29)-(31)给出。

围岩中的应力场通过叠加三个不同的分量来表述：开挖前的初始应力、开挖后的应力状态以及支护安装后的应力。完整的表达式由方程(32)-(34)给出。

围岩和衬砌的切向位移类似地获得，它们的表达式由方程(35)和(36)给出。

在[0°–90°]区间内分析了三种接触模式下的接触应力分布，利用了圆形隧洞的固有对称性。为保持一致性并便于分析应力分布模式，图中的应力以应力p的形式呈现。计算采用侧压力系数λ = 0.5和粘聚力c = 0.15 MPa。如图2所示，当摩擦系数f_r = 0且粘聚力c = 0时，界面接触应力与光滑接触解完全对应。相反，当f_r ≥ 0.59时，结果完全收敛到完全粘结接触解。对于中间的f_r值，解始终位于这两种极限情况之间。

图3的进一步分析表明，衬砌与围岩之间的最大相对滑动发生在45°位置，即拱肩处，滑动幅度随着f_r的增加而逐渐减小。当水平初始应力小于垂直分量时，45°拱肩处的径向应力σ_r保持不变。在光滑接触条件下（f_r = 0），在拱顶（90°)和拱底（0°)观察到峰值径向应力。然而，当摩擦系数超过0.2时，应力峰值向拱腰移动。随着f_r继续增加，径向应力围绕隧洞周边的分布向扁平椭圆模式转变，其特征是拱顶应力减小和拱腰应力升高。

剪应力τ_rθ和切向位移差呈现反比关系：剪应力的增加对应于相对位移差的减小，两个参数都在拱肩附近达到最大值。当摩擦系数达到临界阈值f_rl = 0.59时，条件|τ_rθ ^p| < f_r|σ_r ^p| + c在整个界面得到完全满足，表明完全过渡到粘结接触行为。因此，在了解界面摩擦系数和粘聚力的情况下，工程师可以可靠地选择适当的接触模型进行设计计算。

3 验证与分析

3.1 解析结果的数值验证

通过三种接触模型的应力和位移结果验证了解析解的准确性。为了进一步验证这些结果，采用ANSYS软件进行数值分析。与历史上设计中使用的传统“荷载-结构模型”不同，本研究实现了包含初始地应力的接触模型。采用“反向应力释放法”来模拟延迟支护安装过程。计算开挖后的等效节点荷载，释放荷载的大小模拟支护延迟效应。在弹性状态下，当位移释放系数为η时，剩余的(1-η)部分荷载作用于衬砌，引起岩-衬相互作用。开挖边界处的最终应力通过叠加初始地应力和释放η倍等效荷载引起的应力来确定。

在ANSYS的几种接触模型中，选择了“面面接触”，因为它具有众多优点，例如支持大滑动和摩擦、提供摩擦应力结果以及允许多种建模控制。在计算过程中，采用PLANE42单元模拟围岩和衬砌组件。围岩与衬砌之间的接触通过ANSYS中的“面面接触”方法实现。具体来说，衬砌表面通过网格划分细分为目标单元（TARGE169），这些目标单元与围岩表面上的接触单元（CONTA171）配对。两者共享相同的实常数集以形成接触对。由于研究涉及摩擦滑移接触问题，KEYOPT(12)关键字选项设置为2以允许相对滑动，两种材料之间的摩擦系数设置为0.5。在本研究中，使用映射函数将ζ平面映射到z平面。因此，理论解相对于数值软件结果顺时针旋转了90°。

ANSYS中模拟示例的边界条件和荷载施加方案如图4所示。所选模型的外部边界尺寸为50.0 m × 50.0 m，开挖隧洞半径为3.0 m。在位置A处施加水平和垂直方向的点约束，而在点B处仅施加垂直约束。在衬砌内边界施加2.0 MPa的表面荷载，计算得到的等效节点荷载施加到接触面上的节点上。

图5显示了衬砌的切向应力和剪应力分布等值线图。在[0°–90°]范围内，切向应力逐渐减小，而剪应力在45°处达到最大值。图6展示了数值解与解析解之间的比较。由于对称性，比较了[0°–90°]范围内的结果。在岩-衬界面处，两种方法获得的接触应力显示出密切的一致性，最大相对误差约3.84%发生在拱顶。