考虑渗透各向异性的水下非圆形隧道渗流场解析

乔 彤，周 建*，张天骄，蒋熠诚

(1.浙江大学 建筑工程学院，浙江 杭州 310012；2.滨海和城市岩土工程研究中心，浙江 杭州 310058)

对于水下隧道各类问题的研究，隧道渗流场是目前一个重要的方面[1-3]，绝大多数学者常将隧道截面假设为圆形以便求解。Zhang等[4]在无限承压含水层井流模型的基础上，利用镜像法与叠加原理推导出了多水平地层条件下圆形隧道涌水量的等效计算公式，分析了隧道埋深、半径和海水深度对涌水量的影响。Meye等[5]基于保角映射法，假定洞周土为各向同性多孔弹性介质，推导了圆形隧道渗流量和周围水头分布公式，研究了不同防水排水类型下的渗流场。在等水压和等水头两种水力边界条件下，Zou等[6]采用坐标变换法，将直角坐标系中稳定渗流控制方程求解问题转化为双极坐标系中的求解问题，简化了圆形隧道渗流量及衬砌外水压力的求解过程。谢寒松[7]以复变函数的反演变换和地下水力学理论为基础，推导了圆形隧道渗流场解析解，并对隧道与溶洞的间距和隧道半径对隧道压力水头的影响进行了分析。张丙强[8]采用镜像法和渗流场叠加法则，推导了含注浆圈与衬砌作用下的圆形隧道渗流场解析解，且其对于非平行的双孔水下隧道同样适用。郭玉峰等[9]基于质量守恒定律和达西定律，采用保角变换法和Schwartz迭代法推导了水下双线平行圆形隧道稳态渗流场的解析式，讨论了隧道间距、埋深和相对大小对双孔平行圆形隧道渗流场水头分布、渗流量的影响。基于Dupuit假设和积分法，Meng等[10]获得了稳定渗流条件下的各向同性垂直分层潜水层中圆形排水隧道渗流场解析解。Tang等[11]基于保角映射复变函数理论推导了考虑渗流各向异性的水下圆形隧洞渗流场解析解，但在映射过程中地层表面将发生变形，因此只能算是一种近似解。基于各向同性稳定渗流假定，Lin等[12]采用保角映射及最优化技术推导了考虑渗流各向异性的水下圆形隧道渗流场的半解析式，并研究了等水压和等水头渗水模式下的隧道渗流量变化情况，但未考虑到衬砌及注浆圈的影响。朱成伟等[13]以保角映射法为基础，结合地下水连续渗流理论及叠加原理，推导了水下双线平行圆形隧道渗流场的解析解，获得了水头分布及隧道涌水量的解析表达式。赵建平等[14]结合双极坐标及镜像法，推导了水下隧道渗流场解析解，并对支护参数取值问题进行了合理分析。徐长节等[15]在渗流各向同性圆形隧洞渗流场解析解的基础上，通过坐标变换及保角映射法，推导了各向异性渗流下深埋圆形隧洞的渗流场解析解，但无限含水层竖井理论的假设与实际不符。

在隧道建设中，矩形、直墙曲拱形、多圆心形等具有高效空间利用率的隧道越来越常见。现有关于水下隧道渗流场的研究常将隧道截面形状假设为圆形，对于非圆形隧道，常采用等效面积法和等效圆周法将其转化为圆形截面的隧道[16-17]。少数关于非圆形隧道渗流场的研究也未考虑到渗流各向异性的影响，关于联合考虑渗流各向异性及非圆形隧道的研究未见报道。

基于此，在地下水渗流力学的基础上，采用坐标变换、保角映射、最优化技术，求解稳定渗流连续方程，给出了在渗透各向异性半无限多孔介质中水下非圆形隧道的渗流量及衬砌外水头的解析推导式。本文提供的方法适用于考虑洞周土渗透各向异性的矩形、直墙曲拱形、多圆心形等任意非圆形隧道的渗流场解析研究。

1 计算模型

本文研究对象为考虑渗透各向异性的水下非圆形隧道，即：水下隧道衬砌横截面并非圆形，而是类似于马蹄形。地层表面以上水头高度为Hw；衬砌内边界形状控制方程为f1=L(x,y), 衬砌外边界形状控制方程为f2=S(x,y) ，渗透系数为k1；隧道中心到地表距离定义为埋深h；洞周土水平向渗透系数为kx，竖直向渗透系数为ky，与渗透主轴方向一致，且均为常数。基准面建立在地层表面；将地表边界与衬砌外边界之间的区域定义为区域1，将衬砌单独定义为区域2。模型示意图具体如图1所示。

图1 非圆形水下隧道模型示意图Fig.1 Schematic diagram of non-circular underwater tunnel model

图2 坐标转换示意图Fig.2 Schematic diagram of coordinate transformation

图3 区域1保角映射Fig.3 Region 1 conformal mapping

为方便推导，模型做如下基本假设：

1）洞周土体为渗透各向异性介质；

2）隧道处于稳定渗流状态；

3）水流服从达西定律；

4）衬砌与土体接触面为密实接触界面，在纵向上满足平面应变条件；

5）排水通过衬砌实现均匀渗水。

按照等水头渗水模式模拟隧道衬砌的排水方式[18]，假设衬砌内水压力为0，则隧道衬砌内缘水头为：

式中，H为渗流域某点水头，h为隧道埋深。

2 解析解推导过程

由于区域1为渗透各向异性介质，区域2衬砌支护为渗透各向同性介质，问题比较复杂。为方便求解，将区域1与区域2分开计算，再按照不同区域渗流量相等原则联立求解[19]。

2.1 区域1渗流场求解

2.1.1 区域1坐标转换

根据地下水渗流的连续性原理，可求得均质、渗透各向异性含水层的2维稳定渗流基本微分方程为：

通过坐标变换将式（2）转化为Laplace方程形式，以方便利用现存Laplace方程求解方法求解区域1渗流场。本文假设渗透系数kx、ky为常数，不发生变化，式（2）两边同除以ky得：

令

式中，X、Y为转换后G平面的横纵坐标。

可得在G平面上的渗流连续方程：

由于内边界形状不规则，在G平面中直接求解式（5）难以实现，可借助保角变换公式将区域1转换为同心圆环进行求解。

2.1.2 区域1保角映射

区域1主要由两个边界条件控制，可采用如式（6）所示的复变函数式[21]，该式可将区域1映射为同心圆环域，将地表边界映射为半径为1的外圆，内边界=s′(X,Y) 映 射为一半径为 α（α ＜1）的内圆。

式中：z为原平面点； ς为像平面点， ς=ξ+iη=ρeiθ，ξ、η 为复变函数w(ς) 的实部和虚部，i为虚数单位，ρ为极径， θ 为极角； µ 为实常数，k=1,2,3,···,µ；a、βk为待定实常数。当给出具体的地表及内边界控制方程=s′(X,Y)时，也可以确定上述值，具体求解过程如下。

1）在地表边界：将地层表面映射为半径为1的圆。

首先，将坐标改写：

将式（7）代入式（6），分离实部和虚部可得：

通过式（8）可知，坐标原点被映射为 ξ=-1，η=0 ；无穷远点映射为 ξ=1 ， η=0；地层表面其余各点也在单位圆上一一对应。

2）在区域1内边界：将内边界映射为半径为 α的圆。

将式（9）代入式（6），分离实部和虚部可得在内边界：

2.1.3 区域1渗流连续方程求解

映射后的圆环域内外半径分别为α、1，边界条件不变， ς平面满足渗流连续方程：

其通解形式为：

式中，H为水头，C1、C2为待定常数，C3、C4为待定常数族，t为常数。代入边界条件：

式中，h1为衬砌外缘水头。

最终得区域1内的水头表达式：

式中，H1为区域1内任一点水头。

在已知洞周土域水头分布的情况下，可以对其进行积分得到区域1的渗流量Q1：

2.2 区域2渗流场求解

2.2.1 区域2保角映射

区域2为马蹄形非圆截面衬砌（图4），包括内外缘两个边界，控制方程为：内边界f1=L(x,y)，外边界f2=S(x,y)。由于是非圆形截面，直接求解困难，本文采用吕爱钟等[22]提出的复变函数式（17）将区域2映射为同心圆环，其中，将衬砌内缘映射为半径为1的内圆，将外缘映射为半径为rs的外圆。

图4 区域2保角映射Fig.4 Region 2 conformal mapping

保角映射公式一般以劳伦级数展开式表示：

式（17）确定后，区域1在式（17）下可映射为像平面的同心圆环域，Laplace方程在像平面依旧成立，可直接进行求解。

2.2.2 区域2渗流连续方程求解

在 λ平面内，Laplace方程形式为：

代入边界条件：

解得衬砌水头H2和渗流量Q2：

式中，H2、Q2分别为区域2的水头和流量。

2.3 区域1、2联立求解

按照层间渗流量相等原则Q1=Q2， 联立式（16）和（21）解得衬砌外水头h1和 隧道渗流量Q：

2.4 存在注浆圈时稳定渗流方程的求解

当注浆圈存在时，其求解思路与第2.2节区域2渗流场的求解思路一致，也即利用复变函数式（17）将注浆圈域转换为同心圆环域，结合边界条件求解渗流连续方程，再由渗流连续性原则与洞周土域、衬砌域联立，不再重复赘述。

3 算例验证

3.1 模型设置

假设水下隧道断面为椭圆形，衬砌内缘长轴l1=4.8 m ，短轴l2=3.0 m ， 厚 0.2 m，埋深20.0 m，水位高于地表20.0 m。以渗透各向异性比（ky为常量）表示洞周土渗透各向异性，取值1、2、3；n越大，表示土体渗透各向异性越强。此外，设衬砌内缘水压力为0，主要参数见表1。

表1 水下隧道计算参数Tab.1 Calculation parameters of an underwater tunnel

表2 区域1参数计算值Tab.2 Area 1 parameter calculation values

3.2 数值模拟

采用COMSOL数值软件进行对比验证，参数取值与表1一致。为消除边界条件的影响，模型宽度取300 m，高度取120 m；采用稳态渗流模式，流体密度取ρ=1 000 kg/m3，孔隙率取0.3；水力边界条件设模型底部及左右侧为不透水边界，上表面为自由透水边界；此外，将有限元网格划分设为加密模式，保证计算结果的收敛性，划分2 971个有限元单元，模型如图5所示。

3.3 解析公式参数计算

3.3.1 区域1参数计算

由本文水下隧道模型参数可知，区域1的内边界形状为长轴为5.0 m，短轴为3.2 m的椭圆，埋深20.0 m，根据第2.2节可计算a、 α 、 βk。经试算， βk取3项即可获得较高的精确度。具体参数计算值见表2。

3.3.2 区域2参数计算

区域2是外边界长短轴分别为5.0、3.2 m，内边界长短轴分别为4.8、3.0 m的椭圆形衬砌。通过式（17）将其转换为同心圆环。其中，对于一般的隧道截面形状只需要取3～4项即可满足精度要求[24]。取v=4，代入计算程序得到：rs=1.051 3，R=3.999 8，Ck值见表3。

表3 v=4时Ck计算结果Tab.3 Ck calculation results when v=4

最终得到映射函数为：

经过上述映射函数式，可将内、外边界映射为半径为1.000 0、1.051 3 m的圆环，且可保证内、外边界获得较高的精确度，最大误差不超过0.263 8 m。

3.4 验证与分析

图6为隧道渗流量随衬砌渗透系数的变化曲线。由图6可知，取有限项v的解析解会与数值模拟解产生一定的误差，但整体上理论解和数值解变化趋势一致，误差很小，拟合度较高，证明了本文解的正确性。观察图6可知，隧道渗流量Q随衬砌渗透系数的增大而增大，但增大趋势渐渐放缓。此外，考虑洞周土渗透各向异性对Q也会产生较大的影响，且衬砌渗透性越好，这种影响越明显。例如：在k1=1×10-7m/s条件下，当n=1时，解析值Q为8.14 m3·d-1·m-1；当n=2、3时，Q为12.67、16.00 m3·d-1·m-1，相比n=1分别增加 55.65%、96.56%。在k1=5×10-8m/s条件下，当n=1时，解析值Q为6.86 m3·d-1·m-1；n=2、3时，Q为9.81、11.70 m3·d-1·m-1，相比n=1分别增加了43.00%、70.55%。从计算结果看，渗透各向异性对隧道渗流量的影响显而易见。与此同时，当衬砌渗透性很低时，n值对Q的影响大大降低，当k1≤1×10-9m/s时，n值的影响已经很小了。由此可见，是否考虑洞周土渗透各向异性的影响要根据水下隧道洞周土及支护结构的渗透性具体分析，不能一概而论。

图6 隧道渗流量随衬砌渗透系数的变化曲线Fig.6 Variation curves of tunnel seepage flow with lining permeability coefficient

图7为水下隧道渗流场水头整体分布规律和隧道衬砌外水头随衬砌渗透性的变化规律。图7（a）为当n=1，k1= 1×10-7m/s时，水下隧道渗流场水头分布。观察图7（a）发现，本文计算的水头分布规律与数值模拟结果基本符合，验证了本文解析法的正确性。图7（b）为衬砌外水头随衬砌渗透系数变化曲线，纵坐标为衬砌外水头与全水头的比（h1/H全，H全为全水头）。由图7（b）可知，衬砌外水头的解析解与数值解基本吻合，且衬砌外水头随着k1的变化规律与隧道渗流量随k1的变化规律一致。当衬砌渗透性很低时，衬砌外缘承担的水头很高，接近于全水头值，这对衬砌的抗压强度提出了很高的要求，必要时需进行注浆加固以减小衬砌外水头带来的水压力。

图7 本文解析法与数值法水头计算结果对比Fig.7 Comparison of analytical and numerical methods for hydraulic head calculation

渗透各向异性比n也会对衬砌外水头产生一定的影响。在图7（b）中，当k1=5×10-8m/s时，n=1的衬砌外水头为全水头的32%，n=2、3的衬砌外水头为全水头的45%、54%。可见，若不考虑洞周土渗透各向异性的影响有可能低估衬砌外水头值，将给衬砌的强度设计带来隐患。

目前，在隧道工程中对于非圆形隧道常采用等效面积或者效周长法将非圆形隧道转化为圆形隧道，等效圆半径计算公式为[25]：

等效周长法：

等效面积法：

式中，r3、r4为等效周长法和等效面积法得到的隧道等效半径，l、S为隧道周长和面积。

本文隧道的内缘长短轴分别为4.8 m、3.0 m，计算周长为26.05 m，面积为45.24 m2，按照等效周长法和等效面积法计算得到r3=4.15 m、r4=3.80 m。同样地，对于衬砌外缘，长短轴分别为5.0、3.2 m。按照等效周长法和等效面积法计算得到的等效圆半径分别为r5=4.35 m、r6=4.00 m。渗流量计算公式采用文献[25]中提供的公式，具体不再展开。

n=1时，使用本文方法与等效周长法和等效面积法所得隧道渗流量计算值随衬砌渗透系数的变化规律如图8所示。从图8可以看出，等效周长法渗流量解整体上大于本文解，等效面积法渗流量解整体上小于本文解，本文解与数值解更吻合，说明本文方法计算的准确性相比等效法更高。此外，当衬砌渗透性增强时，等效法与本文方法计算的渗流量差异将扩大，说明在衬砌渗透性较好或者无衬砌支护的毛洞隧道工况下采用本文方法求解隧道渗流量更有意义。

图8 本文方法与等效法下渗流量对比Fig.8 Comparison of seepage flow between the proposed method and the equivalent methods

衬砌外水压力的准确计算对于支护结构设计具有重要意义。按照式（22）计算出衬砌外水头值后，可以按式（27）计算出衬砌外水压力值：

式中，p为水压力，H总为总水头，y为位置水头， γw为水的重度。

表4为本文方法与等效法在n=1（渗透各向同性）情况下，拱腰和拱顶点处衬砌外水压力计算值。

表4 本文方法与等效法衬砌外水压力比较Tab.4 Comparison of external water pressure between the proposed method and the equivalent method kРa

从表4可以看出，本文方法和等效法在拱腰点处衬砌水压力差异较小，而在拱顶点差异较大。主要原因是等效法将非圆形隧道等效为圆形隧道，改变了衬砌上各点的位置坐标，进而改变了各点的位置水头，导致衬砌外水压力出现较大差异。

综合来看，本文的计算方法在隧道渗流量、衬砌上各点水压力的计算方面相对于等效法更准确，更符合实际，在隧道排水系统设置、衬砌抗压强度设计方面具有重要工程应用价值。

4 结 论

1）通过坐标变换及两个保角映射复变函数式，推导了考虑洞周土渗透各向异性的水下非圆形隧道渗流场解析解，并通过连续性原则联立求解隧道渗流量及衬砌外水头。本文提供的方法适用于考虑洞周土渗透各向异性的任意形状水下隧道渗流场的解析研究。

2）通过对某水下隧道进行算例验证，利用COMSOL数值模拟软件建立模型，与解析解值进行对比，从隧道渗流量及衬砌外水头两方面验证了本文解的正确性。此外，通过与等效周长法和等效面积法计算的隧道渗流量和衬砌外水压力对比发现，本文方法的计算值更准确，更符合实际情况。

3）洞周土渗透各向异性对水下隧道渗流场有重要影响。通过对隧道渗流量及衬砌外水头的分析可知，不考虑洞周土渗透各向异性将会低估上述两者的值，给隧道工程的设计带来安全隐患，因此有必要将洞周土渗透各向异性的影响纳入隧道排水系统及衬砌强度设计之中。

下一步将分析注浆圈对渗透各向异性下非圆形隧道渗流场的影响。