水工圆形隧洞围岩衬砌摩擦滑动接触的新解法1)

尹崇林 吕爱钟

(华北电力大学水电与岩土工程研究所，北京 102206)

引言

隧洞是广泛应用于水电、交通、矿山及军事工程领域的常用地下结构.近来，随着我国经济的持续发展及高新技术的不断应用，我国的隧洞工程得到了前所未有的迅速发展，支护-围岩相互作用，以及如何保证隧洞围岩的稳定性是我们要研究的基本问题[1-2].

利用解析法求解隧洞工程问题是一种基本的分析方法.而由Muskhelishvili[3]提出的平面弹性复变函数方法是求解隧洞围岩和衬砌中应力以及位移的常用解析法，通过这种方法容易获得无支护隧洞的应力和位移解[4-5].而在实际工程中，常在隧洞内设置衬砌支护以保证隧洞的安全，因此研究支护-围岩的相互作用更是十分重要，何川等[6]以及吴顺川等[7]分别利用D-P准则和莫尔库伦屈服准则，获得了围岩和衬砌相互作用的弹塑性解析解.

Hasanpour 等[8]、Ramoni 等[9]以及Son[10]都对支护和岩土材料之间的相互作用问题进行了研究.在求解有支护的深埋隧洞问题的解时，可以将其简化为平面应变无限域问题，于学馥等[11]、Wang 等[12]以及Lu 等[13-14]已经获得了各向同性弹性岩体中圆形隧洞衬砌和围岩相互作用的解析解.Lu等[15-16]考虑了更复杂的孔形，以直墙半圆拱形隧洞为例获得了非圆形隧洞的解析解.Bobet[17]和王少杰等[18]考虑了更复杂的岩体材料，分别视岩体为横观各向同性和正交各向异性的弹性体，获得了非圆形隧洞的应力和位移的解析解.

隧洞在施加衬砌支护之后，围岩和衬砌接触并相互作用，在以往的解析研究中，大都考虑了完全接触和光滑接触两种理想的接触方式.文献[12-13,15,17-18]都假定围岩与衬砌之间的接触为完全接触，即认为接触界面上法向的径向应力和径向位移及环向的剪切应力和切向位移都是连续的.这种接触假定围岩与衬砌之间非常粗糙，接触面可以承受很大的摩擦力，不允许围岩与衬砌之间产生相对滑动；而文献[11,14,16]在求解过程中将围岩与衬砌之间的接触视为光滑接触，即认为接触界面上法向的径向应力和径向位移仍然是连续的，即围岩和衬砌在法线方向不能相互脱开，但假定接触面充分光滑，在接触面的切线方向不能承受摩擦力，即认为环向方向的剪切应力为零，这种接触允许围岩与衬砌之间产生相对滑动.Atkinson 和Eftaxiopoulos[19]在求解套管井和胶结井水力压裂问题的解析解时，在水泥和岩石的接触面上考虑了这种光滑接触条件.高永涛等[20]考虑双层厚壁圆筒之间的光滑接触并获得了非均布荷载作用下内外壁的应力解析解.储昭飞等[21]在求解非静水应力场中圆形隧道黏弹性解析解时，其中也考虑了衬砌和围岩之间的光滑接触.

完全接触和光滑接触是人为假定的两种理想情形，实际上两种材料的接触很可能并非为这两种理想情形，例如，在建立连接结构接触界面的非线性力学模型时，王东等[22]就将名义的光滑平面视作凹凸不平的粗糙面，考虑了微凸体的黏滑摩擦行为.围岩和衬砌接触时，它们之间并非完全光滑，也并非可以承受任意大的摩擦力，当围岩与衬砌之间的剪应力大于所能承受的最大静摩擦力后，将发生切向滑动.库仑摩擦模型作为刚体的摩擦模型，因其简单和实用而被广泛应用于工程分析中.该模型认为切向摩擦力的数值不能超过法向压力和摩擦系数的乘积.Meschke 等[23]建立计算模型以模拟桩土相互作用时隧道的机械化开挖，其中就利用了库仑摩擦模型来模拟桩-土的相互作用.Cavalieri 等[24-25]利用增广拉格朗日法求解涉及摩擦的接触问题时，也利用库仑摩擦定律建立了计算模型.苏宗贤等[26]和杨钊等[27]利用库仑摩擦模型分别模拟了管片和围岩以及复合衬砌中内外层衬砌之间的接触.吕爱钟等[28]在库仑摩擦模型的基础上利用复变函数方法和最优化方法获得了圆形隧洞围岩和衬砌摩擦滑动接触的解析解，并与有限元软件ANSYS 所得结果进行了对比，验证了其方法合理性.

本文同样将通过库仑摩擦模型模拟接触面上的切向力学行为，但不同于文献[28]，本文在优化过程中减少了设计变量的个数，极大地简化了优化模型.在弹性接触分析中，库仑摩擦模型可表述如下

式中，τrθ接触面上的剪应力，σr为接触面法向压力，fr为接触界面综合摩擦系数.该表达实质上是静摩擦力的屈服条件或接触面的滑移条件，将切向力与径向压力联系起来.这种接触条件最符合实际情况，即为摩擦滑动接触.本文基于平面应变假定，将围岩和衬砌视为各向同性线弹性体，考虑支护滞后效应，利用复变函数方法推导出深埋圆形隧洞在原始地应力和衬砌内部静水压力共同作用下围岩与衬砌在摩擦滑动接触情况下的解析解，为隧洞围岩和衬砌接触的支护设计计算提供理论基础.

1 基本原理及方法

如图1所示，为在原始地应力和静水压力共同作用下圆形衬砌支护隧洞.在开挖之前，围岩无穷远处水平和竖直方向的原始地应力分别为：，λ 为侧压力系数，定义拉应力为负，压应力为正.隧洞衬砌的内外半径分别为R0和R1，定义径向和x方向的夹角为θ，以逆时针为正.

图1 原始地应力和静水压力共同作用下深埋圆形衬砌隧洞Fig.1 Lined circular tunnel under in situ stresses and uniform hydrostatic pressure

1.1 摩擦滑动接触的准则及数学模型

当围岩和衬砌间的接触为摩擦滑动接触时，由于围岩和衬砌之间的摩擦作用，接触面可以传递剪力，但允许围岩与衬砌之间的切向相对滑动，本文将通过最优化方法来处理这样的问题.当围岩与衬砌之间的剪应力大于所能承受的最大静摩擦力后，围岩和衬砌间接触面间将产生相对滑动，接触面两侧，即围岩一侧和衬砌一侧的切向位移将产生间断，即在r=R1的接触面上，会存在uθ1uθ2的区域.本文作者提出在保证围岩与衬砌之间剪应力小于或者等于接触面所能提供的摩擦力前提下，将接触面间产生最小滑移量的状态作为衬砌的真实工作状态.这样的准则可以用最优化的数学模型表示为

其中，F(X)=|uθ1−uθ2|，bj(X)0 为库仑摩擦模型表示的不等式组.

1.2 仅开挖作用时引起的围岩位移

如图2所示在隧洞开挖之后，假设没有支护，围岩会完成相应的位移.此时围岩内任一点的径向位移ur和切向位移uθ可以表示为[13]

图2 开挖引起的围岩位移Fig.2 Surrounding rock displacement caused by excavation

式中，R1为无支护隧洞的半径，z是一个复数且z=R1eiθ，κ1=3-4µ1，G1=E1/[2(1+µ1)]，µ1，E1和G1分别是岩石的泊松比、弹性模量和剪切模量.式(3)即隧洞开挖后无支护时围岩的全部位移.

1.3 求解围岩和衬砌解析函数的基本方程

对于圆形隧洞，其应力分量和位移分量极易求得[3]

其中，ϕ0(z)和ψ0(z)是关于围岩或者衬砌的两个解析函数.若能求出满足边界条件的围岩和衬砌的两个解析函数，则应力分量和位移分量均可以利用上述结果得出.

隧洞开挖后进行支护，围岩和衬砌相互作用.围岩对应的两个解析函数为[13]

同时，衬砌对应的两个解析函数可以用洛朗级数表示为

式(9)和式(10)中的系数ck，dk，ek，fk，gk和hk都是待求的实常数.

通过衬砌内边界L1的应力边界条件和围岩衬砌接触面L2的接触条件可以建立求解解析函数系数的基本方程.

用σr2，σθ2，τrθ2分别表示衬砌的径向应力、环向应力和剪应力.在衬砌内边界L1上，作用有大小为p0的静水压力，则其应力边界条件可以表示为(式中z=R0eiθ)

用σr1，σθ1，τrθ1分别表示围岩的径向应力，环向应力和剪应力.在围岩衬砌接触面L2上，围岩和衬砌的径向应力和剪应力连续，其应力连续条件可以表示为

本文考虑支护滞后的隧洞开挖过程，η 为位移释放系数，如果η=0，即隧洞在开挖后立即进行支护，围岩在支护之前没有发生位移.明显这不符合工程实际，即使在开挖后立即安装支护，也不可避免地在围岩中发生了一部分变形.假设隧洞位移在完成了最大位移ur+iuθ的η(0η1)倍后，再进行支护，此时围岩完成的位移为η(ur+iuθ).设置衬砌以后，围岩和衬砌相互作用，支护限制了部分围岩位移的产生，此时围岩的位移表示为ur1+iuθ1.而在围岩作用下，衬砌的位移可以表示为ur2+iuθ2.

在围岩衬砌接触面L2上，根据法向位移连续条件可得，其中z=R1eiθ

式中κ2=3-4µ2，G2=E2/[2(1+µ2)].µ2，E2和G2分别是衬砌的泊松比、弹性模量和剪切模量.

将相应的z值代入式(9)∼式(12)，式(15)∼式(17)，利用幂级数解法，比较θ的同幂次项系数，可得到关于系数ck，dk+2，ek，fk+2，gk+2和hk的方程，其中k1.联立所有方程，整理后发现当k2时，ck，dk+2，ek，fk+2，gk+2和hk都为0，只有c1，d1，d3，e1，f1，f3，g1，g3和h1为非零实数，可得方程(18)∼式(25)

这样利用边界条件就得到了式(18)∼式(25)共8 个线性方程.8 个方程，9 个未知量c1，d1，d3，e1，f1，f3，g1，g3，h1，只通过这8 个方程，是不可能求解9 个未知量的，如何求解这样的问题是本文的关键所在.

本文将用围岩一侧和衬砌一侧的切向位移间断值uθ1−uθ2的大小来衡量接触面间的滑移量，在式(2)中用目标函数F(X)=|uθ1−uθ2| 表示.由式(8)可得(式中z=R1eiθ)

由式(28)和式(29)可得

由式(30)可以看出，对任意θ，当|{}|内的值达到最小值时，所产生的位移间断值都最小，只是对于不同θ，位移间断值的最小值不同而已.所以取目标函数F为

使式(31)达到最小的一组c1，d1，d3，e1，f1，f3，g1，g3和h1，可以保证接触面间产生最小的滑移量.在接触面，必须满足库仑摩擦模型，式中τrθ和σr分别为接触面上的剪应力和径向正应力，由于在接触面上应力的连续性，因此τrθ，σr既可取τrθ1，σr1，也可取τrθ2，σr2.本文取前者，则有

对任意的θ，式(32)都应该成立.由于问题的对称性，只考虑θ∈[0,π]的情形即可.将[0,π]分为m等分(本文取180)，则式(32)可化归为m+1 个约束条件

2 摩擦滑动接触的解法

2.1 最优化问题的求解

文献[28]使用混合罚函数方法来进行优化计算[29-30]，其数学模型为

按照文献[28]的求解过程，其将式(18)∼式(25)组成的8 个方程作为优化模型中的等式约束函数aj(X)，而将式(33)组成的不等式组作为不等式约束条件，然后将等式约束函数aj(X)和不等式约束函数bj(X)分别构成惩罚项加到目标函数F(X)构成一个新的无约束的称为罚函数的目标函数p(x,r)，其中r称为罚因子.这样，就把原目标函数F(X)的约束极小值问题转为求罚函数p(x,r)的无约束极小值问题.其罚函数的具体形式为

式中下标集合I1，I2定义为

其中X(0)是给出的原问题的计算初始点.在该模型中，X=[c1,d1,d3,e1,f1,f3,g1,g3,h1]，即所有的变量都需要参与到优化的过程中，计算初始点X(0)=我们知道，设计变量的增加会优化的数学模型更加复杂，从而增加优化的难度，有时甚至得不到满足约束的最优解.如果可以减少设计变量的个数，这将大大地简化优化的数学模型，从而使求解变得容易.

本文同样使用混合罚函数法来进行优化计算，不同的是本文将设计变量的个数从9 个变为1 个，极大地简化了优化模型.具体来说，本文可以将c1，d1，d3，e1，f1，f3，g1，g3和h1其中的任意一个作为设计变量，以h1为例，赋予h1初值，根据式(18)∼式(25)组成方程组即可以求解出c1，d1，d3，e1，f1，f3，g1，g3这8 个变量，再将c1，d1，d3，e1，f1，f3，g1，g3和h1代入如下的优化模型中，即可获得优化结果

式中，bj(X)0 表示的是不等式约束条件，它由式(33)构成，而F(X)即为式(31).在该优化模型中，其罚函数的具体形式简化为

可以看出，较之文献[28]中的罚函数的具体形式，本文的罚函数形式没有二次损失项中的并且计算初始点X(0)仅有，这极大地简化了优化的数学模型.模型中只有h1为设计变量，改变h1的值，使式(36)达到最小值的即为该优化问题的解，此时其他8个变量可以由解方程组得出且分别表示为为了比较两种模型的精确程度，我们以摩擦滑动接触为例，利用新旧两种模型分别计算出不同摩擦系数fr对应的目标函数的值，其对比如表1所示.

表1 改变摩擦系数时两种优化模型所得目标函数值比较Table 1 Objective function values obtained by two optimization models when changing friction coefficient fr

从表1 可以看出，比较不同摩擦系数fr计算出的目标函数值，本文提出的新优化模型所得结果更小，这是因为本文通过解方程的方法可以使式(34)中的aj(X)=0 精确满足，较之文献[28]中的优化模型，该方法更加精确.

同时可以通过访问目标函数的次数来对比两种方法的迭代速度(见表2)，可以看出新优化模型较之旧模型有很大地提升.

表2 两种优化模型访问目标函数的次数Table 2 Times to access F(X)for two optimization models

利用新优化模型求解出的这样一组X=[c1,d1,d3,e1,f1,f3,g1,g3,h1]即为本文问题的解，从而完成了围岩支护后，围岩和衬砌中复势函数的求解.从以上的分析和求解过程可以看出，这种摩擦滑动接触的解法具有一般性，它包含了完全接触和光滑接触两种极限情形.

当fr=0 时，由式(32)表示的不等式约束条件将变为

若令τrθ1=0，则可以获得与式(38)完全相同的方程，即式(38)表示的就是光滑接触的情形.

当fr的值大于某个值时，式(32)表示的不等式约束条件总会满足，式(36)达到最小值(其值为零)时的解为

则式(39)所表示的实际就是uθ1=uθ2，这由式(26)、式(27)可以清晰地看出.由后面的算例可知，满足完全接触条件所需要的fr值大小与侧压力系数λ 密切相关，λ 越接近于1，所需要的fr值越小，理论上当λ=1 时，fr=0 就可以满足完全接触条件.

2.2 衬砌和围岩中的应力

给定R0，R1，p，λ，p0，G1，µ1，G2，µ2，η，可以求出待定的c1，d1，d3，e1，f1，f3，g1，g3，h1.

衬砌中的应力可表示为[13]

围岩中的应力应由开挖前，开挖后，支护后三部分的应力进行迭加.迭加后的应力仍用σr1，σθ1，τrθ1符号表示，可求得[13]

3 分析和讨论

3.1 围岩和衬砌接触面上的接触方式

取计算参数为：p=10.0 MPa，p0=2.0 MPa，R1=3.0 m，R0=2.7 m，µ1=0.25，µ2=0.20，E2/E1=10.0，η=0.20.对于不同的侧压力系数λ以及摩擦系数fr，获得的解析函数的系数也不同.由于问题的对称性，只需要分析θ∈[0◦,90◦]内的结果.

由前面的推导我们知道，本文作者提出的摩擦滑动接触可以计算光滑接触和完全接触两种极限的接触工况，并且这与摩擦系数fr的取值相关.取侧压力系数λ=0.5，分析摩擦系数fr取不同值时接触面上接触应力σr，τrθ和切向位移间断值|uθ1−uθ2|的变化规律，并将其与完全接触[13]和光滑接触[14]的结果进行对比.

从图3 可以看出，接触面上一旦发生滑动，除了对称点0◦和90◦，其余各点都会发生相对滑动，且在45◦处有最大值.当摩擦系数fr为零时，接触上发生最大的相对滑动，当fr0.64 时，接触面上的切向位移间断值|uθ1−uθ2|基本为零，即uθ1≈uθ2，并且随着fr的增大，|uθ1−uθ2|逐渐减小.

图3 接触面上切向位移间断值|uθ1−uθ2|的分布Fig.3 Distributions of|uθ1−uθ2|on the interface

如图4(a)所示，接触面上径向应力σr在45◦处为一定值，而当fr=0 即光滑接触时，σr的绝对值分别在0◦和90◦处取得极小值和极大值.当fr较大时，σr的绝对值在[0◦,90◦]范围内单调递减.图4(b)表明随着fr的增大，接触面上的剪应力τrθ的绝对值也随之增大，且都在45◦处取得最大值.当fr0.64时，图4 中的σr，τrθ不再变化，且与完全接触的结果一致，结合图3，fr0.64 时，uθ1≈uθ2，此时得到的解即为完全接触的解.

图4 接触面上接触应力随摩擦系数fr 变化的分布规律Fig.4 Distributions of the contact stresses on the interface for different fr

我们已经确定利用摩擦滑动接触解法可以求解完全接触和光滑接触两种极限接触问题，而且通过推导已知当fr的值大于某个值时，总能满足完全接触的条件，称这个值为阈值确定阈值对判断接触面上的接触方式很有帮助.文献[28]通过反复试算的方式讨论了不同弹性模量和侧压力系数所确定的阈值，但是这种方法不够精确，本文尝试利用一种精确的方法来确定阈值

表3 改变位移释放系数η 时不同侧压力系数λ所确定的阈值Table 3 Values of for various η and λ

表3 改变位移释放系数η 时不同侧压力系数λ所确定的阈值Table 3 Values of for various η and λ

表4 改变衬砌厚度R1−R0时不同侧压力系数λ所确定的阈值Table 4 Values of for various R1−R0and λ

表4 改变衬砌厚度R1−R0时不同侧压力系数λ所确定的阈值Table 4 Values of for various R1−R0and λ

表5 改变水压力p0时不同侧压力系数λ 所确定的阈值Table 5 Values of for various p0and λ

表5 改变水压力p0时不同侧压力系数λ 所确定的阈值Table 5 Values of for various p0and λ

同时，为了同文献[28]中所得结果进行比较，取R0=2.5 m，侧压力系数分别为0.2，0.5，0.8，1.0，在不同弹模比值下所得阈值如表6 所示.

表6 改变弹模比值E2/E1时不同侧压力系数λ所确定的阈值Table 6 Values of for various E2/E1and λ

表6 改变弹模比值E2/E1时不同侧压力系数λ所确定的阈值Table 6 Values of for various E2/E1and λ

3.2 衬砌和围岩各边界上切向应力的变化规律

在隧洞的支护设计中，各边界上的切向应力变化规律十分重要，所以在这一节，我们取不同的摩擦系数fr，分别分析它们的变化对衬砌内外边界的切向应力以及围岩开挖边界上切向应力(σθ1)的影响.3.1 节已经分析过当λ=0.5，满足完全接触条件的阈值=0.64，因此，在本节的讨论中，摩擦系数fr的最大取值为0.64，其余参数不变化.

由图5 可以看出，当λ=0.5，不管摩擦系数如何改变，所有的切向应力都为压应力.由图5(a)和图5(b)可以看出，随着fr的增大，衬砌内边界的切向应力最大压应力增大，最小压应力减小，在[0◦,90◦]内，压应力随着角度的增大而减小，其变化范围增大.在衬砌外边界上，随着fr的增大，切向应力变化范围减小.只是当fr较小且趋于0时，从0◦到90◦，压应力随着角度的增大而增大，分别在0◦和90◦取得最小和最大压应力；当fr较大时，结论恰好相反.由图5(c)可得，保持λ=0.5 不变，在[0◦,90◦]内，随着角度的增大，σθ1为压应力单调减小.随着fr的增大，σθ1在90◦出现的最小压应力增大，在0◦出现的最大压应力减小.图5 表明摩擦滑动接触条件下所得衬砌内外边界上的切向应力以及围岩开挖边界上的切向应力σθ1都介于光滑接触和完全接触的结果之间.

图5 各边界上切向应力随摩擦系数fr 变化的分布规律Fig.5 Distributions of the tangential stresses on the 3 boders for different fr

4 结论

本文提出了更加符合实际情况的摩擦滑动接触条件来描述实际工程中围岩和衬砌之间的接触问题.以库仑摩擦模型模拟围岩衬砌之间的接触，当接触面上发生相对滑动之后，将接触面上产生最小滑移量的状态视为衬砌的真实工作状态，并以此为基础，利用最优化方法建立了摩擦滑动接触情形下新的解法.在利用混合罚函数法求解解析函数的系数时，减少了设计变量的个数，去掉了罚函数二次损失项中的极大地简化了最优化解法的数学模型.

同时提出了能够精确获得满足完全接触条件的摩擦系数阈值的公式，并且得到了不同位移释放系数，不同衬砌厚度，不同水压力大小以及不同弹模比值情况下改变侧压力系数时所获得的阈值为判断围岩和衬砌接触方式提供了理论基础.摩擦滑动接触条件下所得的解都介于完全接触和光滑接触之间，并且能同时可以求解光滑接触和完全接触两种极限情况，具有一般性.特别地，在衬砌内边界上，摩擦滑动接触所得最大切向应力介于两种极限接触之间，且光滑接触所得最大切向应力最小，这说明在实际工程中，尽量减小衬砌和围岩之间的摩擦，可以增加衬砌的承载能力.